Secular orbital dynamics of the innermost exoplanet of the $\upsilon$-Andromedæ system

Abstract

We introduce a quasi-periodic restricted Hamiltonian to describe the secular motion of a small-mass planet in a multi-planetary system. In particular, we refer to the motion of $\upsilon$-And $\textbf{b}$ which is the innermost planet among those discovered in the extrasolar system orbiting around the $\upsilon$-Andromedæ A star. We preassign the orbits of the Super-Jupiter exoplanets $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$ in a stable configuration. The Fourier decompositions of their secular motions are reconstructed by using the well-known technique of the (so-called) frequency analysis and are injected in the equations describing the orbital dynamics of $\upsilon$-And $\textbf{b}$ under the gravitational effects exerted by those two external exoplanets (that are expected to be major ones in such an extrasolar system). Therefore, we end up with a Hamiltonian model having $2+3/2$ degrees of freedom; its validity is confirmed by the comparison with several numerical integrations of the complete four-body problem. Furthermore, the model is enriched by taking into account also the effects due to the relativistic corrections on the secular motion of the innermost exoplanet. We focus on the problem of the stability of $\upsilon$-And $\textbf{b}$ as a function of the parameters that mostly impact on its orbit, that are the initial values of its inclination and the longitude of its node (as they are measured with respect to the plane of the sky). In particular, we study the evolution of its eccentricity, which is crucial to exclude orbital configurations with high probability of (quasi)collision with the central star in the long-time evolution of the system. Moreover, we also introduce a normal form approach, that is based on the complete average of our restricted model with respect to the angles describing the secular motions of the major exoplanets. Therefore, our Hamiltonian model is further reduced to a system with 2 degrees of freedom, which is integrable because it admits a constant of motion that is related to the total angular momentum. This allows us to very quickly preselect the domains of stability for $\upsilon$-And $\textbf{b}$, with respect to the set of the initial orbital configurations that are compatible with the observations.

1 Introduction

The $\upsilon$-Andromedæ system was the first ever to be discovered among the ones that host at least two exoplanets. In fact, a few years after the discovery of the first exoplanet, the evidence for multiple companions of the $\upsilon$-Andromedæ system was announced (see Mayor and Queloz 1995; Butler et al. 1999, respectively). In particular, the observations made with the detection technique of the radial velocity (hereafter, RV) revealed the existence of orbital objects with three different periods, 4.6, 241 and 1267 days, which revolve around $\upsilon$-Andromedæ A, that is the brightest star of a binary hosting also the red dwarf $\upsilon$-Andromedæ B. Such exoplanets were named $\upsilon$-And $\textbf{b}$, $\textbf{c}\,$ and $\textbf{d}\,$ in increasing order with respect to their distance from the main star. Since $\upsilon$-Andromedæ B is very far with respect to these other bodies (i.e., $\sim 750\,$AU), it is usual to not consider its negligible gravitational effects when the dynamical behavior of the planetary system orbiting around $\upsilon$-Andromedæ A is studied.

None of the present detection methods allows us to know all the orbital parameters of an extrasolar planet. For instance, the RV technique does not provide any information about both the inclination and the longitude of node. In the $\upsilon$-Andromedæ case, these two orbital elements were measured (although with rather remarkable error bars) for both $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$ thanks to observational data taken from the Hubble Space Telescope (see McArthur et al. 2010). The information provided by such an application of astrometry significantly complemented the knowledge about this extrasolar system; in fact, it has led to the evaluation of the masses of $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$ (ranging in $13.98^{+2.3}_{-5.3}\>M_J$ and $10.25^{+0.7}_{-3.3}\>M_J$, respectively) and of their mutual inclination ($29.9^\circ \pm 1^\circ$). It is well known that only minimum limits for the masses can be deduced by observations made with the RV method (due to the intrinsic limitations of such a technique). Moreover, the mutual inclination between planetary orbits plays a crucial role for what concerns the stability of extrasolar systems (see, e.g., Volpi et al. 2019). Thus, the orbital configuration of $\upsilon$-Andromedæ is probably one of the most accurately known among the extrasolar multi-planetary systems which have been discovered so far.

The question of the orbital stability of $\upsilon$-Andromedæ is quite challenging. Numerical integrations revealed that unstable orbits are frequent. Moreover, these extensive explorations allowed to locate four different regions of initial values of the orbital parameters (consistent with all observational constraints) yielding dynamically stable orbital configurations for the three planets of the system (see Deitrick et al. 2015). All these sets of parameters correspond to values of the mass of $\upsilon$-And $\textbf{c}$ that are relatively small, in the sense that they are much closer to the lower bound of the range $13.98^{+2.3}_{-5.3}\>M_J$ than to the upper one. On the other hand, according to a numerical criterion inspired from normal form theory and introduced in Locatelli et al. (2022), the most robust orbital configurations correspond to the largest possible value of the mass of $\upsilon$-And $\textbf{c}$ in the above range. In the vicinity of the initial conditions giving rise to the most robust orbital configuration, the existence of KAM tori for the dynamics of a secular three-body problem including $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$ was proved in a rigorous computer–assisted way (see Caracciolo et al. 2022).

In the present paper, we aim to extend the study of the stability to the orbital dynamics of $\upsilon$-And $\textbf{b}$, still adopting a hierarchical approach. In the case of the particular extrasolar system under consideration, this means that we assume the mass of $\upsilon$-And $\textbf{b}$ so small¹ with respect to the ones of $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$, that the motion of the innermost planet $\textbf{b}$ can be modeled with a good approximation via a restricted four-body problem. More precisely, in order to study the dynamical behavior, we preassign the secular motions of the Super-Jupiter exoplanets $\textbf{c}$ and $\textbf{d}$ in correspondence to the quasi-periodic orbit that is expected to be the most robust. The Fourier decompositions of the secular motions of $\textbf{c}$ and $\textbf{d}$ are reconstructed by using the well-known technique of the (so-called) frequency analysis (see, e.g., Laskar 2005) and are injected in the equations describing the orbital dynamics of $\upsilon$-And $\textbf{b}$, under the gravitational effects exerted by those two external exoplanets. This way to introduce a quasi-periodic restricted model has been recently used to study the long-term dynamics of our Solar System (see Mogavero and Laskar 2022; Hoang et al. 2022).

The advantage of introducing a secular quasi-periodic restricted Hamiltonian looks evident. In our present case (referring to the $\upsilon$-Andromedæ system with planets $\textbf{b}\,$, $\textbf{c}\,$, $\textbf{d}$), we start with a Hamiltonian model having 9 degrees of freedom, ending up with a simpler one with $2+3/2$ degrees of freedom, where the short periods are dropped. This explains why numerical explorations of the restricted Hamiltonian model are much faster. Our main purpose is further promoting this procedure, in such a way to introduce a new simplified model with just two degrees of freedom. Indeed, in the present paper, we show that this can be done by adopting a suitable normal form approach, which is so accurate to produce an integrable Hamiltonian that can be used to efficiently characterize the stability domain with respect to the unknown orbital parameters of $\upsilon$-And $\textbf{b}$, i.e., the inclination and the longitude of the node.

The present work is organized as follows. In Sect. 2, frequency analysis is used to reconstruct the Fourier decompositions of the secular motions of the outer exoplanets $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$. In Sect. 3, the secular quasi-periodic restricted Hamiltonian model (with $2+3/2$ degrees of freedom) is introduced and validated through the comparison with several numerical integrations of the complete four-body problem, hosting planets $\textbf{b}\,$, $\textbf{c}\,$, $\textbf{d}\,$ of the $\upsilon$-Andromedæ system. The double normalization procedure allowing to perform a sort of averaging which further simplifies the model is described with a rather general approach in Sect. 4. In Sect. 5, this normal form procedure is applied to the quasi-periodic restricted Hamiltonian, in such a way to derive an integrable model with 2 degrees of freedom describing the secular orbital dynamics of $\upsilon$-And $\textbf{b}\,$. Such a simplified model is used to study $\upsilon$-And $\textbf{b}$ stability domain in the parameters space of the initial values of the inclination and the longitude of node. All this computational procedure is repeated in Sect. 6 starting from a version of the secular quasi-periodic restricted Hamiltonian model which includes also relativistic corrections; this allows us to appreciate the effects on the orbital dynamics due to general relativity.

2 Determination of the outer planets motion

To prescribe the orbits of the giant planets $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$, we start from the Hamiltonian of the three-body problem (hereafter, 3BP) in Poincaré heliocentric canonical variables, using the formulation based on the reduced masses $\beta _2\,$, $\beta _3\,$, that is

$$\begin{aligned} \mathcal {H}=&\sum _{j=2}^{3}\left( \frac{{{\varvec{p}}}_j\cdot {{\varvec{p}}}_j}{2\, \beta _j} -\frac{\mathcal {G}\, m_0\, m_j}{r_j}\right) +\frac{{{\varvec{p}}_2}\cdot {{\varvec{p}}_3}}{ m_0} -\frac{\mathcal {G}\, m_2\, m_3}{|{{\varvec{r}}_2}-{{\varvec{r}}_3}|}\, , \end{aligned}$$

(1)

where $m_0$ is the mass of the star, $m_{j}\,$, ${{\varvec{r}}}_j\,$, ${{\varvec{p}}}_j\,$, $j=2,3$, are the masses, astrocentric position vectors and conjugated momenta of the planets, respectively, $\mathcal {G}$ is the gravitational constant and $\beta _j=m_0 m_j/(m_0+m_j)$, $j=2,3$, are the reduced masses. Let us remark that, in the following, we use the indexes 2 and 3, respectively, for the inner ($\upsilon$-And $\textbf{c}$) and outer ($\upsilon$-And $\textbf{d}$) planets between the giant ones, while the index 1 is used to refer to $\upsilon$-And $\textbf{b}$

In order to set up a quasi-periodic restricted model for the description of the motion of $\upsilon$-And $\textbf{b}$, we need to characterize the motion of the giant planets; this can be done through the frequency analysis method, starting from the numerical integration of the complete 3BP Hamiltonian, reported in Eq. (1) (i.e., before any expansions and averaging²). Thus, we numerically integrate the complete Hamiltonian (1) using a symplectic method of type $\mathcal {SBAB}_3\,$, which is described in Laskar and Robutel (2001). As initial orbital parameters for the outer planets, we adopt those reported in Table 1, corresponding to the most robust planetary orbit compatible with the observed data available for $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$ (see McArthur et al. 2010), according to the criterion of “minimal area” explained in Locatelli et al. (2022).

Table 1 Chosen values of the masses and of the initial orbital parameters for $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$, compatible with the observed data available, as reported in McArthur et al. (2010)

Having fixed as initial orbital parameters the ones described in Table 1, it is possible to compute their correspondent values in the Laplace reference frame (i.e., the invariant reference frame orthogonal to the total angular momentum vector ${{\varvec{r}}}_2\times {{\varvec{p}}}_2+{{\varvec{r}}}_3\times {{\varvec{p}}}_3\,$) and to perform the numerical integration of the full 3BP corresponding to these initial values. Then, it is possible to express the discrete results produced by the numerical integrations in the canonical Poincaré variables $(\xi _j, \eta _j)$, $(P_j, Q_j)$ (momenta coordinates) given by³

$$\begin{aligned} \begin{aligned}&\xi _j=\sqrt{2\Gamma _j}\cos (\gamma _j)=\sqrt{2\Lambda _j}\,\sqrt{1-\sqrt{1-\textrm{e}_j^2}}\cos (\varpi _j)\,, \\&\eta _j=\sqrt{2\Gamma _j}\sin (\gamma _j)=-\sqrt{2\Lambda _j}\,\sqrt{1-\sqrt{1-\textrm{e}_j^2}}\sin (\varpi _j)\,,\qquad j=1,\,2,\,3\,,\quad \\&P_j=\sqrt{2\Theta _j}\cos (\theta _j)=2\sqrt{\Lambda _j}\,\root 4 \of { 1-\textrm{e}_j^2}\,\sin \left( \frac{\textrm{i}_j}{2}\right) \cos ( \Omega _j)\,,\\&Q_j=\sqrt{2\Theta _j}\sin (\theta _j)=-2\sqrt{\Lambda _j}\,\root 4 \of { 1-\textrm{e}_j^2}\,\sin \left( \frac{\textrm{i}_j}{2}\right) \sin ( \Omega _j)\, \end{aligned} \end{aligned}$$

(2)

where $\Lambda _j=\beta _j\sqrt{\mu _j a_j}\,$, $\beta _j=m_0 m_j/(m_0+m_j)\,$, $\mu _j=\mathcal {G}\left( m_0+m_{j}\right) \,$, and $\textrm{e}_j\,$, $\textrm{i}_j\,$, $\omega _j\,$, $\Omega _j$, $\varpi _j=\omega _j+\Omega _j$ refer, respectively, to the eccentricity, inclination, argument of the periastron, longitudes of the node and of the periastron of the j-th planet.

However, the numerical integrations do not allow to obtain a complete knowledge of the motion laws $t\mapsto (\xi _j(t), \eta _j(t))$, $t\mapsto (P_j(t), Q_j(t))$ ($j=2,3$), producing only discrete time series made by sets of finite points computed on a regular grid in the interval $[0,T]\,$. The computational method of frequency analysis (hereafter, FA) allows, however, to reconstruct with a good accuracy the motion laws by using suitable continuous in the time variable t functions. This has been done recently in Mogavero and Laskar (2022) and Hoang et al. (2022), in order express the motion of the Jovian planets of our Solar System as a Fourier decomposition including just a few of the main terms. In the present Section, we basically follow that approach; therefore, here we limit ourselves to report some definitions which are essential in order to make our computational procedure well definite (see, e.g., Laskar 2005 for an introduction and a complete exposition about FA). We consider analytic quasi-periodic motion laws $t\mapsto z(t)$. This means that the function $z:\mathbb {R}\mapsto \mathbb {C}$ admits the following Fourier series decomposition:

$$\begin{aligned} z(t)=\sum _{{{\varvec{k}}}\in \mathbb {Z}^n}a_{{{\varvec{k}}}}e^{i({{\varvec{k}}}\cdot {{\varvec{\omega }}}t + \vartheta _{{\varvec{k}}})}, \end{aligned}$$

(3)

where ${{\varvec{\omega }}}\in \mathbb {R}^n$ is the so-called fundamental angular velocity vector, while $a_{{\varvec{k}}}\in \mathbb {R}_+\cup \lbrace 0 \rbrace$ and $\vartheta _{{\varvec{k}}}\in \mathbb {T}\,$, $\forall \>{{\varvec{k}}}\in \mathbb {Z}^n$; moreover, the sequence $\{a_k\}_{{{\varvec{k}}}\in \mathbb {Z}^n}$ is assumed to satisfy an exponential decay law, i.e., $|a_{{{\varvec{k}}}}|\le c\, e^{-|{{\varvec{k}}}|\sigma }$ $\forall \,{{\varvec{k}}}\in \mathbb {Z}^n$ with c and $\sigma$ real positive parameters. The FA allows us to find an approximation of z(t) of the following form:

$$\begin{aligned} z(t)\simeq \sum _{s=1}^{\mathcal {N}_C}a_{s;T}e^{i(\nu _{T}^{(s)}t+\vartheta _{s;T})} \end{aligned}$$

(4)

where ${\mathcal {N}_C}$ is the number of components we want to consider. Equation (4) is an approximation of the motion z(t) in the sense that if ${\mathcal {N}_C}\rightarrow +\infty$ and $T\rightarrow +\infty \,$, the right-hand side of (4) converges to (3). Moreover, $a_{s;T}\in \mathbb {R_+}\,$ and $\vartheta _{s;T}\in [0,2\pi )$ are called, respectively, the amplitude and the initial phase of the s-th component, while $\nu _{T}^{(s)}$ is a local maximum point of the function

$$\begin{aligned} \nu \mapsto \mathcal {A}(\nu )=\frac{1}{T}\left| \int _{0}^{T} z(t)\,e^{-i\nu t}\,\mathcal {W}(t)\, \textrm{d}t \right| \,, \end{aligned}$$

(5)

where $\mathcal {W}$ is a suitable weight function such that ${\int _{0}^{T}\mathcal {W}(t)\,\textrm{d}t=T}\,$. Following Laskar (2005), we use the so-called Hanning filter, defined (in [0, T]) as $\mathcal {W}(t)=1-\cos [\pi (2t/T-1)]\,$.

The numerical integration of the 3BP (Eq. 1) produces a discretization of the signals⁴$t\mapsto \xi _j(t)+i\eta _j(t)$ and $t\mapsto P_j(t)+i Q_j(t)\,$, that are $\lbrace (\xi _j(s\Delta ), \eta _j(s\Delta ))\rbrace _{s=0}^{\mathcal {N}_P}$ and $\lbrace (P_j(s\Delta ), Q_j(s\Delta ))\rbrace _{s=0}^{\mathcal {N}_P}\,$ ($j=2,3$), where $s = 0, \ldots ,\mathcal {N}_P$ and the sampling time is $\Delta = T /\mathcal {N}_P\,$. These discretizations allow to compute (by numerical quadrature) the integral in (5) and, consequently, a few of local maximum points of the function (5) considering, as z(t), $\xi _2(t)+i\eta _2(t)\,$, $\xi _3(t)+i\eta _3(t)\,$, $P_2(t)+iQ_2(t)$ and $P_3(t)+iQ_3(t)\,$.

Then, we use the FA to compute a quasi-periodic approximation of the secular dynamics of the giant planets $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$, i.e.,

$$\begin{aligned} \begin{aligned}&\xi _j(t)+i\eta _j(t)\simeq \sum _{s=1}^{\mathcal {N}_C}A_{j,s}e^{i({{\varvec{k}}}_{j,s}\cdot {{\varvec{\theta }}}(t)+\vartheta _{j,s})}\,, \\&P_j(t)+iQ_j(t)\simeq \sum _{s=1}^{\widetilde{\mathcal {N}}_C}\tilde{A}_{j,s}e^{i(\tilde{{{\varvec{k}}}}_{j,s}\cdot {{\varvec{\theta }}}(t)+\tilde{\vartheta }_{j,s})}\,, \end{aligned} \end{aligned}$$

(6)

$j=2,3\,$, where the angular vector

$$\begin{aligned} {{\varvec{\theta }}}(t)=(\theta _3(t), \theta _4(t), \theta _5(t))=(\omega _3\,t, \omega _4\, t, \omega _5\, t):={{\varvec{\omega }}}\,t \end{aligned}$$

(7)

depends linearly on time and ${{\varvec{\omega }}}\in \mathbb {R}^3$ is the fundamental angular velocity vector whose components are listed in the following:

$$\begin{aligned} \begin{aligned}&\omega _3=-2.4369935819462266 \times 10^{-3}\, , \\&\omega _4=-1.04278712796661375 \times 10^{-3}\, , \\&\omega _5=4.88477275490260560 \times 10^{-3}\, . \end{aligned} \end{aligned}$$

(8)

Hereafter, the secular motion of the outer planets $t\mapsto (\xi _j(t), \eta _j(t), P_j(t), Q_j(t))$, $j=2,3$, is approximated as it is written in both the r.h.s. of formula (6). The numerical values of the coefficients which appear in the quasi-periodic decompositions⁵ of the motions laws are reported in Tables 2, 3, 4 and 5.

In order to verify that the numerical solutions are well approximated by the quasi-periodic decompositions computed above, we compare the time evolution of the variables $\xi _2$, $\xi _3$, $\eta _2$, $\eta _3$, $P_2$, $P_3$, $Q_2$, $Q_3$ as obtained by the numerical integration and by the FA. Figure 1 shows that the quasi-periodic approximations nearly perfectly superpose to the plots of the numerical solutions.

Table 2 Decomposition of the signal $\xi _2(t)+i\,\eta _2(t)\,$ as it is provided by the FA

Table 3 Decomposition of the signal $\xi _3(t)+i\,\eta _3(t)\,$ as it is provided by the FA

Table 4 Decomposition of the signal $P_2(t)+i\,Q_2(t)\,$ as it is provided by the FA

Table 5 Decomposition of the signal $P_3(t)+i\,Q_3(t)\,$ as it is provided by the FA

3 The secular quasi-periodic restricted (SQPR) Hamiltonian

Having preassigned the motion of the two outer planets $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}\,$, it is now possible to properly define the secular model for a quasi-periodic restricted four-body problem (hereafter, 4BP). We start from the Hamiltonian of the 4BP, given by

Fig. 1

Dynamical behavior of the variables $\xi _2$, $\xi _3$, $\eta _2$, $\eta _3$, $P_2$, $P_3$, $Q_2$, $Q_3$ (their definition is reported in (2)) as it is computed by numerical integration of the complete (non secular) three-body problem and by the quasi-periodic approximation provided by the FA; the corresponding plots are in blue and red, respectively

$$\begin{aligned} \mathcal {H}_{4BP}= \sum _{j=1}^{3}\left( \frac{{{\varvec{p}}}_j\cdot {{\varvec{p}}}_j}{2\, \beta _j} -\frac{\mathcal {G}\, m_0\, m_j}{r_j}\right) +\sum _{1\le i< j \le 3}\frac{{{\varvec{p}}_i}\cdot {{\varvec{p}}_j}}{ m_0} -\sum _{1\le i < j \le 3}\frac{\mathcal {G}\, m_i\, m_j}{|{{\varvec{r}}_i}-{{\varvec{r}}_j}|}\,. \end{aligned}$$

(9)

We recall that the so-called secular model of order one in the masses is given by averaging with respect to the mean motion angles, i.e.,

$$\begin{aligned} \mathcal {H}_{sec}({{\varvec{\xi }}}, {{\varvec{\eta }}}, {{\varvec{P}}}, {{\varvec{Q}}})= \>\int _{\mathbb {T}^3}\!\!\!\!\frac{\mathcal {H}_{4BP}({{\varvec{\Lambda }}}, {{\varvec{\lambda }}}, {{\varvec{\xi }}}, {{\varvec{\eta }}}, {{\varvec{P}}}, {{\varvec{Q}}})}{8\pi ^3}\, \textrm{d} \lambda _1 \textrm{d}\lambda _2 \textrm{d}\lambda _3\,. \end{aligned}$$

(10)

In the l.h.s. of the equation above, we disregard the dependence on the actions ${{\varvec{\Lambda }}}$, because in the secular approximation of order one in the masses their values $\Lambda _j=\beta _j\sqrt{\mu _j\,a_j}\,$, $j=1,2,3$, are constant. Due to the d’Alembert rules (see, e.g., Murray and Dermott 1999; Morbidelli 2002), it is well known that the secular Hamiltonian can be expanded in the following way:

$$\begin{aligned} \mathcal {H}_{sec}({{\varvec{\xi }}}, {{\varvec{\eta }}}, {{\varvec{P}}}, {{\varvec{Q}}})= \sum _{s=0}^{\mathcal {N}/2}\sum _{\begin{array}{c} |{{\varvec{i}}}|+|{{\varvec{l}}}|+\quad \,\,\,\\ |{{\varvec{m}}}|+|{{\varvec{n}}}|=2s \end{array}}\!\!\!\!\!\! c_{{{\varvec{i}}},{{\varvec{l}}},{{\varvec{m}}},{{\varvec{n}}}} \prod _{j=1}^{3}\xi _j^{i_j}\eta _j^{l_j} P_{j}^{m_j} Q_j^{n_j}\,, \end{aligned}$$

(11)

where $\mathcal {N}$ is the order of truncation in powers of eccentricity and inclination. We fix $\mathcal {N}=8$ in all our computations.

Since we aim at describing the dynamical secular evolution of the innermost planet $\upsilon$-And $\textbf{b}\,$, it is sufficient to consider the interactions between the two pairs $\upsilon$-And $\textbf{b}$, $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{b}$, $\upsilon$-And $\textbf{d}$. In more details, let $\mathcal {H}_{sec}^{\mathfrak {i} -\mathfrak {j}}$ be the secular Hamiltonian derived from the three-body problem for the planets $\mathfrak {i}$ and $\mathfrak {j}$, averaging with respect to the mean longitudes $\lambda _\mathfrak {i}\,$, $\lambda _\mathfrak {j}\,$. Its expansion writes as

$$\begin{aligned} \mathcal {H}^{\mathfrak {i}-\mathfrak {j}}_{sec}(\xi _{\mathfrak {i}}, \eta _{\mathfrak {i}}, P_{\mathfrak {i}}, Q_{\mathfrak {i}},\xi _{\mathfrak {j}}, \eta _{\mathfrak {j}}, P_{\mathfrak {j}}, Q_{\mathfrak {j}})=\sum _{s=0}^{\mathcal {N}/2}\sum _{\begin{array}{c} |{{\varvec{i}}}|+|{{\varvec{l}}}|+\quad \,\,\,\\ |{{\varvec{m}}}|+|{{\varvec{n}}}|=2s \end{array}}\!\!\!\!\!\! c_{{{\varvec{i}}},{{\varvec{l}}},{{\varvec{m}}},{{\varvec{n}}}}\prod _{j=\mathfrak {i},\mathfrak {j}}\xi _j^{i_j}\eta _j^{l_j} P_{j}^{m_j} Q_j^{n_j}\,. \end{aligned}$$

(12)

Therefore, a restricted non-autonomous model which approximates the secular dynamics of $\upsilon$-And $\textbf{b}$ can be defined by considering the terms

$$\begin{aligned}{} & {} \mathcal {H}^{1-2}_{sec}(\xi _1, \eta _1, P_1, Q_1,\xi _2(t), \eta _2(t), P_2(t), Q_2(t)) \\{} & {} \quad + \mathcal {H}^{1-3}_{sec}(\xi _1, \eta _1, P_1, Q_1,\xi _3(t), \eta _3(t), P_3(t), Q_3(t)) \,, \end{aligned}$$

where $\xi _2(t)$, $\eta _2(t)$, ...$P_3(t)$, $Q_3(t)$ are replaced with the corresponding quasi-periodic approximations written in both the r.h.s. appearing in formula (6). Let us stress that, having prescribed the motion of the two outermost planets $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}\,$, at this stage the Hamiltonian $\mathcal {H}_{sec}^{2-3}$ does not need to be reconsidered; indeed, it would introduce additional terms that disappear in the equations of motion (see formula (15)).

We can finally introduce the quasi-periodic restricted Hamiltonian model for the secular dynamics of $\upsilon$-And $\textbf{b}\,$; it is given by the following $2+3/2$ degrees of freedom Hamiltonian:

$$\begin{aligned}{} & {} \mathcal {H}_{sec,\, 2+\frac{3}{2}}({{\varvec{p}}},{{\varvec{q}}},\xi _1, \eta _1, P_1, Q_1 )= \omega _3\,p_3 + \omega _4\,p_4 + \omega _5\,p_5 \nonumber \\{} & {} \quad +\mathcal {H}^{1-2}_{sec}(\xi _1, \eta _1, P_1, Q_1, \xi _2({{\varvec{q}}}), \eta _2({{\varvec{q}}}), P_2({{\varvec{q}}}), Q_2({{\varvec{q}}})) \nonumber \\{} & {} \quad +\mathcal {H}^{1-3}_{sec}(\xi _1, \eta _1, P_1, Q_1,\xi _3({{\varvec{q}}}), \eta _3({{\varvec{q}}}), P_3({{\varvec{q}}}), Q_3({{\varvec{q}}}))\,, \end{aligned}$$

(13)

where the pairs of canonical coordinates referring to the planets $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$ (that are $\xi _2\,$, $\eta _2\,$, ...$P_3\,$, $Q_3$) are replaced by the corresponding finite Fourier decomposition written in formula (6) as a function of the angles ${{\varvec{\theta }}}$, renamed⁶ as ${{\varvec{q}}}$, i.e.,

$$\begin{aligned} {{\varvec{q}}}=(q_3,\,q_4,\,q_5):= (\theta _3,\,\theta _4,\,\theta _5)={{\varvec{\theta }}}\,. \end{aligned}$$

(14)

Let us focus on the summands appearing in the first row of (13), i.e., the Hamiltonian term ${{\varvec{\omega }}}\cdot {{\varvec{p}}}\,$, where ${{\varvec{\omega }}}$ is the fundamental angular velocity vector (defined in formula (8)) and ${{\varvec{p}}}=(p_3, p_4, p_5)$ is made by three so-called dummy variables, which are conjugated to the angles ${{\varvec{q}}}$. The role they play is made clear by the equations of motion for the innermost planet, which write in the following way in the framework of the restricted quasi-periodic secular approximation:

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q_3}=\partial \mathcal {H}_{sec,\,2+\frac{3}{2}}/\partial p_3=\omega _3\\ \dot{q_4}=\partial \mathcal {H}_{sec,\,2+\frac{3}{2}}/\partial p_4=\omega _4\\ \dot{q_5}=\partial \mathcal {H}_{sec,\,2+\frac{3}{2}}/\partial p_5=\omega _5\\ \dot{\xi }_{1}=-\partial \mathcal {H}_{sec,\,2+\frac{3}{2}} /\partial \eta _1= -\partial \big ( \mathcal {H}^{1-2}_{sec} + \mathcal {H}^{1-3}_{sec} \big ) /\partial \eta _1 \\ \dot{\eta }_{1}=\partial \mathcal {H}_{sec,\,2+\frac{3}{2}} /\partial \xi _1= \partial \big ( \mathcal {H}^{1-2}_{sec} + \mathcal {H}^{1-3}_{sec} \big ) /\partial \xi _1 \\ \dot{P}_{1}=-\partial \mathcal {H}_{sec,\,2+\frac{3}{2}} /\partial Q_1= -\partial \big ( \mathcal {H}^{1-2}_{sec} + \mathcal {H}^{1-3}_{sec} \big ) /\partial Q_1 \\ \dot{Q}_{1}=\partial \mathcal {H}_{sec,\,2+\frac{3}{2}}/\partial P_1= \partial \big ( \mathcal {H}^{1-2}_{sec} + \mathcal {H}^{1-3}_{sec} \big ) /\partial P_1 \end{array}\right. }\,. \end{aligned}$$

(15)

Due to the occurrence of the term ${{\varvec{\omega }}}\cdot {{\varvec{p}}}$ in the Hamiltonian $\mathcal {H}_{sec,\,2+\frac{3}{2}}\,$, the first three equations admit ${{\varvec{q}}}(t)={{\varvec{\omega }}}t$ as a solution, in agreement with formulæ (7) and (14). This allows to reinject the wanted quasi-periodic time dependence in the Fourier approximations $\xi _2({{\varvec{q}}})$, $\eta _2({{\varvec{q}}})$, $\ldots$ $P_3({{\varvec{q}}})$, $Q_3({{\varvec{q}}})$. As a matter of fact, we do not need to compute the evolution of $({p}_3(t),{p}_4(t),{p}_5(t))$ because they do not exert any influence on the motion of $\upsilon$-And $\textbf{b}\,$; they are needed just if one is interested in checking that the energy is preserved, because it is given by the evaluation of $\mathcal {H}_{sec,\,2+\frac{3}{2}}\,$.

We also recall that, in order to produce a restricted quasi-periodic secular model, it is possible to apply the closed-form averaging, which is compared in Mastroianni (2023) with the computational method that is adopted here and is based on the expansions in Laplace coefficients. Finally, we emphasize what is discussed below.

Remark 3.1

The Hamiltonian $\mathcal {H}_{sec,\,2+\frac{3}{2}}$ is invariant with respect to a particular class of rotations. Thus, it admits a constant of motion that could be reduced, so to decrease⁷ by one the number of degrees of freedom of the model.

In order to clarify the statement above, it is convenient to introduce a complete set of action-angle variables, defining two new pairs of canonical coordinates $\xi _1=\sqrt{2\Gamma _1}\cos (\varpi _1)$, $\eta _1=\sqrt{2\Gamma _1}\sin (-\varpi _1)$, $P_1=\sqrt{2\Theta _1}\cos (\Omega _1)$, $Q_1=\sqrt{2\Theta _1}\sin (-\Omega _1)$, referring to a pair of orbital angles of $\upsilon$-And $\textbf{b}$, i.e., $\varpi _1$ and $\Omega _1\,$, that are the longitudes of the pericenter and of the node, respectively (see definition (2) of the Poincaré canonical variables). Thus, it is possible to verify the following invariance law:

$$\begin{aligned}{} & {} -\frac{\partial \,\mathcal {H}_{sec,\,2+\frac{3}{2}}}{\partial \xi _1}\, \,\frac{\partial \,\xi _1}{\partial \varpi _1}- \frac{\partial \,\mathcal {H}_{sec,\,2+\frac{3}{2}}}{\partial \eta _1}\, \frac{\partial \,\eta _1}{\partial \varpi _1}- \frac{\partial \,\mathcal {H}_{sec,\,2+\frac{3}{2}}}{\partial P_1}\, \frac{\partial P_1}{\partial \Omega _1}- \frac{\partial \,\mathcal {H}_{sec,\,2+\frac{3}{2}}}{\partial Q_1}\, \frac{\partial Q_1}{\partial \Omega _1} \\{} & {} \quad + \frac{\partial \,\mathcal {H}_{sec,\,2+\frac{3}{2}}}{\partial q_3}+ \frac{\partial \,\mathcal {H}_{sec,\,2+\frac{3}{2}}}{\partial q_4}+ \frac{\partial \,\mathcal {H}_{sec,\,2+\frac{3}{2}}}{\partial q_5}\,=\,0\,. \end{aligned}$$

Therefore, $p_3+p_4+p_5+\Gamma _1+\Theta _1$ is preserved; such a quantity, apart from an inessential additional constant, is equivalent to the total angular momentum.

The above invariance law is better understood, recalling that $q_3$ and $q_4$ correspond to the longitudes of the pericenters of $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$, respectively, while $q_5$ refers to the longitude of the nodes of $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$ (that in the Laplace frame, determined by taking into account only these two exoplanets are opposite one to the other). This identification is due to the way we have determined $(q_3,q_4,q_5)$ by decomposing some specific signals of the secular dynamics of the outer exoplanets (this is made by using the frequency analysis, as it is explained in Sect. 2). Thus, the aforementioned invariance law is due to the fact that the dynamics of the model we are studying does depend just on the pericenters arguments of the three exoplanets and on the difference between the longitude of the nodes of $\upsilon$-And $\textbf{b}$ and $\upsilon$-And $\textbf{c}$, i.e., $\Omega _1-\Omega _2=\Omega _1-\Omega _3-\pi$. Therefore, the Hamiltonian is invariant with respect to any rotation of the same angle that is applied to all the longitudes of the nodes; as it is well known, by Noether theorem, this is equivalent to the preservation of the total angular momentum.

Table 6 Available data for the orbital parameters of the exoplanet $\upsilon$-And $\textbf{b}$. The values above are reported from Table 13 of McArthur et al. (2010)

3.1 Numerical validation of the SQPR model

In order to validate our secular quasi-periodic restricted (hereafter SQPR) model describing the dynamics of $\upsilon$-And $\textbf{b}$, we want to compare the numerical integrations of the complete 4BP with the ones of the equations of motion (15). Let us recall that the chosen values of parameters and initial conditions for the two outer planets are given in Table 1. For what concerns the orbital elements of the innermost planet $\upsilon$-And $\textbf{b}$, both the inclination $\textrm{i}_1$ and the longitude of the ascending node $\Omega _1$ are unknown (see, e.g., McArthur et al. 2010). The available data for $\upsilon$-And $\textbf{b}$ are reported in Table 6. Among the possible values of the initial orbital parameters of $\upsilon$-And $\textbf{b}$, we have chosen $a_1\,$, $\textrm{e}_1$, $M_1$ and $\omega _1$ as in the stable prograde trial PRO2 described in Deitrick et al. (2015). They are reported in Table 7 and are compatible with the available ranges of values appearing in Table 6. Let us recall that the dynamical evolution of the SQPR model does not depend on the mass of $\upsilon$-And $\textbf{b}$; therefore, the choice about its value is not reported in Table 7. For what concerns the unknown initial values of the inclination and of the longitude of nodes, we have decided to vary them so as to cover a 2D regular grid of values $(\textrm{i}_1(0),\,\Omega _1(0))\in I_\textrm{i}\times I_\Omega =[6.865^\circ , 34^\circ ]\times [0^\circ ,360^\circ ]$, dividing $I_\textrm{i}$ and $I_\Omega$, respectively, in 20 and 60 sub-intervals; this means that the widths of the grid steps are equal to $1.35675^\circ$ and $6^\circ$ in inclination and longitude of nodes, respectively. Let us recall that the lowest possible value of the interval $I_\textrm{i}$, i.e., $\textrm{i}_2(0)=6.865^\circ$, corresponds to the inclination of $\upsilon$-And $\textbf{c}$. Considering values smaller than $\textrm{i}_2(0)$ could be incoherent with the assumptions leading to the SPQR model we have just introduced; indeed, the factor $1/\sin (\textrm{i}_1(0))$ increases the mass of the exoplanet by one order of magnitude with respect to the minimal one. Therefore, for small values of $\textrm{i}_1(0)$ the mass of $\upsilon$-And $\textbf{b}$ could become so large that the effects exerted by its gravitational attraction on the outer exoplanets could not be neglected anymore. On the other hand, it will be shown in the sequel that the stability region for the orbital motion of $\upsilon$-And $\textbf{b}$ nearly completely disappears for values of $\textrm{i}_1(0)$ larger than $34^\circ$. These are the reasons behind our choice about the lower and upper limits of $I_\textrm{i}$.

We emphasize that the study of the stability domain, as it is deduced by the numerical integrations, can help us to obtain information about the possible ranges of the unknown values $(\textrm{i}_1(0),\Omega _1(0))$. Moreover, the comparisons between the numerical integrations of the complete 4BP and the ones of the SQPR model aim to demonstrate that the agreement is sufficiently good so that it becomes possible to directly work with the latter Hamiltonian model, that has to be considered easier than the former, because the degrees of freedom are $2+3/2$ instead of 9.

Table 7 Values of the initial orbital parameters for $\upsilon$-And $\textbf{b}$ as they have been selected in the stable prograde trial PRO2 of Deitrick et al. (2015) (Table 3)

3.1.1 Numerical integration of the complete four-body problem

For each pair of values $(\textrm{i}_1(0)\,,\Omega _1(0))\in I_\textrm{i}\times I_\Omega$ ranging in the $20\times 60$ regular grid we have previously prescribed, we first compute the corresponding initial orbital elements of the three exoplanets in the Laplace reference frame, then we perform the numerical integration of the complete 4BP Hamiltonian (9) by using the symplectic method of type $\mathcal {SBAB}_3\,$. Contrary to the SPQR model, the numerical integrations of the 4BP are affected by the mass of $\upsilon$-And $\textbf{b}$; its value is simply fixed in such a way that $m_1=0.674/\sin (\textrm{i}_1(0))\,$.

Fig. 2

Color-grid plots of the maximal value reached by the eccentricity of $\upsilon$-And $\textbf{b}$ (on the left) and by the mutual inclination between $\upsilon$-And $\textbf{b}$ and $\upsilon$-And $\textbf{c}$ (on the right). The maxima are computed during the symplectic numerical integrations of the 4BP which cover a timespan of $10^5$ yr

Fig. 3

Color-grid plots of the maximal value reached by the eccentricity of $\upsilon$-And $\textbf{b}$. The maxima are computed during the symplectic numerical integrations of the 4BP which cover a timespan of $10^5$ yr. The results are obtained by numerical integrations which refer to sets of initial conditions that differ (passing from one panel to another) just because of the choice of the initial values of the mean anomaly; from left to right, the plots refer to $M_1(0)$ equal to $51.4286^\circ$, $205.714^\circ$, and $308.571^\circ$, respectively

The largest value reached by the eccentricity can be considered as a very simple numerical indicator about the stability of the orbital configurations. The maximum eccentricity obtained along each of the $21\times 60$ numerical integrations is reported in the left panel of Fig. 2. In particular, since we are interested in initial conditions leading to regular behavior, i.e., avoiding quasi-collisions, every time that the eccentricity $\textrm{e}_1$ exceeds a threshold value (fixed to be equal to 0.85), in the color-grid plots its maximal value is arbitrarily set equal to one. Moreover, since we expect that $\upsilon$-And $\textbf{c}$ is the most massive planet in that extrasolar system and being it the closest one to $\upsilon$-And $\textbf{b}$, it is natural to focus the attention also on the mutual inclination between $\upsilon$-And $\textbf{b}$ and $\upsilon$-And $\textbf{c}\,$. Let us recall that it is defined in such a way that

$$\begin{aligned} \cos (\textrm{i}_{{mut}_{\textbf{b}\textbf{c}}})= \cos (\textrm{i}_1)\cos (\textrm{i}_2)+\sin (\textrm{i}_1)\sin (\textrm{i}_2)\cos (\Omega _1-\Omega _2)\,; \end{aligned}$$

(16)

therefore, for each numerical integration it is also possible to compute the maximal value reached by $\textrm{i}_{{mut}_{\textbf{b}\textbf{c}}}\,$. The results are reported in the right panel of Fig. 2. In both those panels, the color-grid plots are provided as functions of the initial values of the longitude of nodes $\Omega _1$ and the inclination $\textrm{i}_1\,$, which are reported on the x and y axes, respectively. By comparing the two plots in Figure panels 2a and b, one can easily appreciate that the regions which have to be considered as dynamically unstable, because the eccentricity of $\textrm{e}_1$ can grow to large values, correspond also to large mutual inclinations of the planetary orbits of $\upsilon$-And $\textbf{b}$ and $\upsilon$-And $\textbf{c}$.

We remark that the value of the initial mean anomaly $M_1(0)$ is missing among the available observational data reported in Table 6. As a matter of fact, mean anomalies of exoplanets are in general so poorly known that usually their values are not reported in the public databases.⁸ However, in order to understand if (and up to what extent) the initial value $M_1(0)$ can affect the dynamics of $\upsilon$-And $\textbf{b}\,$, we repeat all the numerical integrations of the 4BP for four different initial values of $M_1\,$, chosen so as to have one of them belonging to each of the quadrants $[0^\circ , 90^\circ ]\,$, $[90^\circ , 180^\circ ]\,$, $[180^\circ , 270^\circ ]\,$ and $[270^\circ , 360^\circ ]\,$. In Fig. 3 we report three examples; in particular, they show the color-grid plots of the maximal value reached by the eccentricity $\textrm{e}_1\,$, taking $M_1(0)$ as $360^\circ /7= 51.4286^\circ$, $4\cdot 360^\circ /7= 205.714^\circ$ and $6\cdot 360^\circ /7=308.571^\circ$, respectively. For what concerns the region $[90^\circ , 180^\circ ]$, let us recall that Fig. 2a refers to $M_1(0)=103.53^\circ$. The comparison between Figs. 2a and 3 shows that the choice of the value of $M_1(0)$ does not seem to produce any remarkable impact on the global structure of the dynamical stability of these exoplanetary orbits.

Moreover, the same conclusion applies also to the increasing factor $1/\sin (\textrm{i}_1(0))$ (with $\textrm{i}_1(0)\in I_\textrm{i}$) which multiplies the minimal mass of $\upsilon$-And $\textbf{b}\,$ in such a way to determine the value of $m_1\,$. In fact, substantial differences are not observed between Figs. 2 and 4.

Fig. 4

Color-grid plots of the maximal value reached by the eccentricity of $\upsilon$-And $\textbf{b}$ (on the left) and by the mutual inclination between $\upsilon$-And $\textbf{b}$ and $\upsilon$-And $\textbf{c}$ (on the right). The maxima are computed during the symplectic numerical integrations of the 4BP which cover a timespan of $10^5$ yr. As the only difference with respect to the numerical integration whose results are reported in Fig. 2, here the mass of $\upsilon$-And $\textbf{b}$ is always kept fixed so as to be equal to its minimal value $m_1=0.674\,$

3.1.2 Numerical integration of the secular quasi-periodic restricted model

We want now to compare the previous results with those found in the SQPR approximation of the four-body problem, performing numerical integrations of the system of equations (15). In order to make these comparisons coherent, also here we consider the data listed in Table 7 as initial conditions for the orbital elements of $\upsilon$-And $\textbf{b}$ which are completed with the values of $(\textrm{i}_1(0)\,,\Omega _1(0))$ ranging in the $20\times 60$ regular grid that covers $I_\textrm{i}\times I_\Omega =[6.865^\circ , 34^\circ ]\times [0^\circ ,360^\circ ]$. At the beginning of the computational procedure, the initial values of the orbital elements are determined in the Laplace reference frame, which is fixed by taking into account only the two outermost planets (i.e., the total angular momentum of the system is given only by the sum of the angular momentum of $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}\,$). Of course, this is made in agreement with our choice to consider a restricted framework, because we are assuming that the mass of $\upsilon$-And $\textbf{b}$ is so small that can be neglected.

For each numerical integration we compute the maximal value reached by the eccentricity $\textrm{e}_1$ and the mutual inclination $\textrm{i}_{{mut}_{\textbf{b}\textbf{c}}}\,$. The results are reported in the color-grid plots of the left and right panels of Fig. 5, respectively. Once again, they are provided as functions of $\Omega _1(0)$ and $\textrm{i}_1(0)\,$, whose values appear on the x and y axes, respectively.

Fig. 5

Color-grid plots of the maximal value reached by the eccentricity of $\upsilon$-And $\textbf{b}$ (on the left) and by the mutual inclination between $\upsilon$-And $\textbf{b}$ and $\upsilon$-And $\textbf{c}$ (on the right). The maxima are computed along the RK4 numerical integrations of the equations of motion (15) of the SQPR model, covering a timespan of $10^5$ yr

Fig. 6

Behavior of the fundamental angular velocity $\nu$ as obtained by applying the frequency map analysis method to the signal $\xi _1(t)+i\eta _1(t)\,$, computed through the RK4 numerical integration of the SQPR model (15), covering a timespan of $1.31072\cdot 10^{5}$ yr. We take, as initial conditions, $(\textrm{i}_1(0),\Omega _1(0))\in \lbrace 6.865^\circ \rbrace \times I_\Omega$ for the left panel and $(\textrm{i}_1(0),\Omega _1(0))\in \lbrace 10.9353^\circ \rbrace \times I_\Omega$ for the right one

Comparing Fig. 2a with  5a and Fig. 2b with 5b, respectively, we can immediately conclude the striking similarity of the color-grid plots, implying the same dependence of the dynamics on the initial values of the orbital elements $\textrm{i}_1(0)$ and $\Omega _1(0)\,$ in either model. In particular, the regions of stability located at the two lateral sides of the plots, where the orbit of $\upsilon$-And $\textbf{b}$ does not become very eccentric, are identical. This occurs also for what concerns the plots of the maximal mutual inclination. However, some discrepancies are evident in the central parts of the panels, i.e., for values of $\Omega _1(0)$ ranging between $90^\circ$ and $270^\circ$. We stress that this lack of agreement between the results provided by the two models is expected in these central regions of the panels. Indeed, let us recall that the SQPR model has been introduced starting from some classical expansions in powers of eccentricities and inclinations. Therefore, it is reasonable to expect a deterioration of the accuracy of the SQPR model in the orbital dynamics depicted in the central regions of the plots where large values of the eccentricity $\textrm{e}_1$ and the mutual inclination are attained. We emphasize that similar remarks about the very strong impact of the initial value $\Omega _1(0)$ on the orbital stability of $\upsilon$-And $\textbf{b}$ can be found in Sect. 4.2 of Piskorz et al. (2017).

A further exploration of the stable and chaotic regions of Fig. 5a can be done by applying the so-called frequency map analysis method (see, e.g., Laskar 1999), in order to study the signal $\xi _1(t)+i\eta _1(t)$ produced by the numerical integration of the system (15), i.e., in the SQPR approximation. We perform the numerical integrations as prescribed at the beginning of the present Section, taking into account only a few values in $I_\textrm{i}\,$ for the initial inclinations, i.e., $\textrm{i}_1(0)=6.865^\circ ,\, 8.22175^\circ ,\,9.5785^\circ ,\,10.9353^\circ$ and $\Omega _1(0)\in I_\Omega \,$. In Fig. 6, we report the behavior of the angular velocity corresponding to the first component of the approximation of $\xi _1(t)+i\eta _1(t)$, as obtained by applying the FA computational algorithm; therefore, this quantity is related to the precession rate of $\varpi _1\,$. As initial value for the inclination $\textrm{i}_1(0)$, we fix $6.865^\circ$ for Fig. 6a and $10.9353^\circ$ for Fig. 6b. We do not report the cases $(\textrm{i}_1(0),\Omega _1(0))\in \lbrace 8.22175^\circ ,\, 9.5785^\circ \rbrace \times I_\Omega \,$, since the behavior of those plots is similar to the ones in Fig. 6.

The situation is well described in Fig. 6b and analogous considerations can be done for Fig. 6a. For what concerns the values of $\Omega _1(0)$ in the range $[0, \sim 50^\circ ]$ and $[\sim 325^\circ , \, 360^\circ ]$ we can observe a regular behavior of the angular velocity $\nu$ which is also monotone with the only exception of the local minimum. According to the interpretation of the frequency map analysis (in light of KAM theory), such a regular regime is due to the presence of many invariant tori which fill the stability region located at the two lateral sides of the plot 5a. Instead for values of $\Omega _1(0)$ in $[\sim 50^\circ , \sim 70^\circ ] \cup [\sim 300^\circ , \, \sim 325^\circ ]$ and $\Omega _1(0)$ in $[\sim 120^\circ , \, \sim 270^\circ ]$ we observe a strongly irregular behavior, which corresponds to the lateral green stripes and the internal green region of Fig. 5a. Thus, they represent chaotic regions in proximity of a secular resonance. Indeed, in Fig. 6b the angular velocity is constant for values of $\Omega _1(0)$ in $[\sim 70^\circ , \, \sim 85^\circ ]$ and $[\sim 280^\circ , \, \sim 300^\circ ]$ (corresponding to part of the blue central stripes of Fig. 5a). More precisely, the value of $\nu$ is equal to $\simeq -1.04\cdot 10^{-3}\,$, that is $\omega _4\,$, i.e., one of the fundamental angular velocities which characterize the quasi-periodic motion of the outer planets (see Eq. (8)). This allows us to conclude that they represent the stable central part of a resonant region.

4 Introduction of a secular model by a normal form approach

This section aims at manipulating the Hamiltonian with normal form algorithms in order to define a new model that is more compact; this allows us to simulate the secular dynamics of $\upsilon$-And $\textbf{b}$ with much faster numerical integrations. In fact, we describe a reduction of the number of degrees of freedom (DOF) of our Hamiltonian model. For such a purpose, we apply two normal form methods: first, we perform the construction of an elliptic torus, hence, we proceed removing the angles ($q_3,q_4,q_5\,$) whose evolution is linearly depending on time. The latter elimination is made by applying a normalization method à la Birkhoff in such a way to introduce a so-called resonant normal form⁹ that includes, at least partially, the long-term effects due to the outer planets motion.

4.1 Normal form algorithm constructing elliptic tori

In Giorgilli et al. (2014), the existence of invariant elliptic tori in 3D planetary problems with n bodies has been proved by using a normal form method which is explicitly constructive. However, such an approach does not look suitable to be directly applied to Hamiltonian secular models, because in this latter case the separation between fast and slow dynamics is lost. Therefore, we follow the explanatory notes (Locatelli et al. 2022), where the algorithm constructing the normal form for elliptic tori is compared with the classical one à la Kolmogorov, which is at the basis of the original proof scheme of the KAM theorem. We first summarize this general procedure leading to the construction of elliptic tori. We then add some comments explaining how this general method can be suitably adapted to our problem.

We start considering a Hamiltonian $\mathcal {H}^{(0)}$ written as follows:

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(0)}&({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})= \mathcal {E}^{(0)}+{{\varvec{\omega }}}^{(0)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(0)}\cdot {{\varvec{I}}}+\sum _{s\ge 0}\sum _{l\ge 3} f_{l}^{(0,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\quad +\sum _{s\ge 1} f_{0}^{(0,s)}({{\varvec{q}}})+\sum _{s\ge 1} f_{1}^{(0,s)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})+\sum _{s\ge 1} f_{2}^{(0,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,, \end{aligned} \end{aligned}$$

(17)

where $\mathcal {E}^{(0)}$ is a constant term, representing the energy, $({{\varvec{p}}},{{\varvec{q}}})\in \mathbb {R}^{n_1}\times \mathbb {T}^{n_1}\,$, $({{\varvec{I}}},{{\varvec{\alpha }}})\in \mathbb {R}^{n_2}_{\ge 0}\times \mathbb {T}^{n_2}$ are action-angle variables and $({{\varvec{\omega }}}^{(0)},{{\varvec{\Omega }}}^{(0)})\in \mathbb {R}^{n_1}\times \mathbb {R}^{n_2}$ is the angular velocity vector. The symbol $f_{l}^{(r,s)}$ is used to denote a function of the variables $({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})\,$, such that l is the total degree in the square root of the actions $({{\varvec{p}}}\,, {{\varvec{I}}})$, s is the index such that the maximum trigonometric degree, in the angles $({{\varvec{q}}},{{\varvec{\alpha }}})\,$, is sK (for a fixed positive integer K) and r refers to the normalization step. In more details, we can say that $f_{l}^{(0,s)}\!\!\in \mathfrak {P}_{l,sK}\,$, which is a class of functions that we introduce as follows.

Definition 4.1

We say that $g\in \mathfrak {P}_{l,sK}$ if $\displaystyle {g\in \!\!\!\!\bigcup _{\begin{array}{c} \widehat{m}\ge 0, \,\widehat{l}\ge 0 \\ 2\widehat{m}+\widehat{l}=l \end{array}}\!\!\!\!\widehat{\mathfrak {P}}_{\widehat{m},\widehat{l},sK}\,,}$ where

$$\begin{aligned}&\widehat{\mathfrak {P}}_{\widehat{m},\widehat{l},sK}=\Bigg \lbrace g:\mathbb {R}^{n_1}\times \mathbb {T}^{n_1}\times \mathbb {R}_{\ge 0}^{n_2}\times \mathbb {T}^{n_2}\rightarrow \mathbb {R}\, :\\&g({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})\!=\!\!\!\!\!\sum _{\begin{array}{c} {{\varvec{m}}}\in \mathbb {N}^{n_1}\\ |{{\varvec{m}}}|=\widehat{m} \end{array}}\sum _{\begin{array}{c} {{\varvec{l}}}\in \mathbb {N}^{n_2}\\ |{{\varvec{l}}}|=\widehat{l} \end{array}}\sum _{\begin{array}{c} {{\varvec{k}}}\in \mathbb {Z}^{n_1}\\ |{{\varvec{k}}}|+|\widehat{{{\varvec{l}}}}|\le sK \end{array}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\sum _{\begin{array}{c} \qquad \quad \widehat{l}_j=-l_j,-l_j+2,\ldots ,l_j\\ j=1,\ldots ,n_2 \end{array}}\!\!\!\!\!\!\!\!\!\!\!c_{{{\varvec{m}}},{{\varvec{l}}},{{\varvec{k}}},\widehat{{{\varvec{l}}}}}\,\,\,{{\varvec{p}}}^{{{\varvec{m}}}}\left( \sqrt{{{\varvec{I}}}}\right) ^{{{\varvec{l}}}}e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}+\widehat{{{\varvec{l}}}}\cdot {{\varvec{\alpha }}}\right) }\Bigg \rbrace \, . \end{aligned}$$

A few remarks about the above definition are in order. First, since we deal with real functions, the complex coefficients must be such that $c_{{{\varvec{m}}},{{\varvec{l}}},-{{\varvec{k}}},-\widehat{{{\varvec{l}}}}}=\bar{c}_{{{\varvec{m}}},{{\varvec{l}}},{{\varvec{k}}},\widehat{{{\varvec{l}}}}}\,$. Moreover, the rules about the integer coefficients vector $\widehat{{{\varvec{l}}}}$ are such that, $\forall \,j=1,\ldots ,n_2\,$, the j-th component of the Fourier harmonic $\widehat{l}_j$ (that refers to the angle $\alpha _j$) must have the same parity with respect to the corresponding degree $l_j$ of $\sqrt{I_j}$ and must satisfy the inequality¹⁰$|\widehat{l}_j|\le l_j\,$.

Let us here emphasize that our SQPR model of the secular dynamics of $\upsilon$-And $\textbf{b}$ can be reformulated in such a way to be described by a Hamiltonian of the type (17); this will be explained in detail at the beginning of Sect. 5.

The following statement plays a substantial role, since it ensures that the structure of the functions $f_{l}^{(r,s)}$ is preserved while the normalization algorithm is iterated.

Lemma 4.1

Let us consider two generic functions $g\in \mathfrak {P}_{l,sK}$ and $h\in \mathfrak {P}_{m,rK}\,$, where K is a fixed positive integer number. Then

$$\begin{aligned} \lbrace g,h \rbrace =L_{h}\,g\,\in \,\mathfrak {P}_{l+m-2,\,(r+s)K}\qquad \forall \,l,\,m,\,r,\,s\,\in \,\mathbb {N}\,, \end{aligned}$$

where we trivially extend the definition 4.1 in such a way that $\mathfrak {P}_{-2,\,sK}=\mathfrak {P}_{-1,\,sK}=\{0\}$ $\forall \ s\in \mathbb {N}$.

The algorithm constructing the normal form for elliptic tori is applied to Hamiltonians of the type (17), where the terms appearing in the second row (namely, $\sum _{s\ge 1}\sum _{l=0}^{2} f_{l}^{(0,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\,$) are considered as the perturbation to remove. Therefore, one can easily realize that such a perturbation must be sufficiently small so that the procedure behaves well as regards convergence. There are general situations where this essential smallness condition is satisfied. For instance, this occurs for Hamiltonian systems in the neighborhood of a stable equilibrium point; in fact, it is possible to prove that, $f_{l}^{(0,s)}=\mathcal {O}(\varepsilon ^s)\,$, where $\varepsilon$ is a small parameter which denotes the first approximation of the distance (expressed in terms of the actions) between the elliptic torus and the stable equilibrium point. The elimination of the small perturbing terms can be done through a sequence of canonical transformations, leading the Hamiltonian in the following final form:

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(\infty )}(\tilde{{{\varvec{p}}}}, \tilde{{{\varvec{q}}}}, \tilde{{{\varvec{I}}}},\tilde{{{\varvec{\alpha }}}})&= \mathcal {E}^{(\infty )}+{{\varvec{\omega }}}^{(\infty )}\cdot \tilde{{{\varvec{p}}}}+{{\varvec{\Omega }}}^{(\infty )}\cdot \tilde{{{\varvec{I}}}}+\sum _{s\ge 0}\sum _{l\ge 3} f_{l}^{(\infty ,s)}(\tilde{{{\varvec{p}}}},\tilde{{{\varvec{q}}}},\tilde{{{\varvec{I}}}},\tilde{{{\varvec{\alpha }}}})\,, \end{aligned} \end{aligned}$$

(18)

with $f_{l}^{(\infty , s)}\in \mathfrak {P}_{l,sK}\,$. Therefore, for any initial conditions of type $({{\varvec{0}}},\tilde{{{\varvec{q}}}}_0, {{\varvec{0}}}, \tilde{{{\varvec{\alpha }}}})$ (where $\tilde{{{\varvec{q}}}_0}\in \mathbb {T}^{n_1}$ and the value of $\tilde{{{\varvec{\alpha }}}}\in \mathbb {T}^{n_2}$ does not play any role¹¹), the motion law $(\tilde{{{\varvec{p}}}}(t), \tilde{{{\varvec{q}}}}(t), \tilde{{{\varvec{I}}}}(t),\tilde{{{\varvec{\alpha }}}}(t))=({{\varvec{0}}},\tilde{{{\varvec{q}}}}_0+{{\varvec{\omega }}}^{(\infty )}t, {{\varvec{0}}}, \tilde{{{\varvec{\alpha }}}})$ is a solution of the Hamilton’s equations related to $\mathcal {H}^{(\infty )}\,$. This quasi-periodic solution (having ${{\varvec{\omega }}}^{(\infty )}$ as constant angular velocity vector) lies on the $n_1$-dimensional invariant torus such that the values of the action coordinates are $\tilde{{{\varvec{p}}}} = {{\varvec{0}}},$ $\tilde{{{\varvec{I}}}}={{\varvec{0}}}\,$.

The generic r-th step of the algorithm is defined as follows. Let us assume that after $r-1$ normalization steps the expansion of the Hamiltonian can be written as

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(r-1)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})&= \mathcal {E}^{(r-1)}+{{\varvec{\omega }}}^{(r-1)}\!\!\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(r-1)}\!\!\cdot {{\varvec{I}}}+\!\sum _{s\ge 0}\sum _{l\ge 3} f_{l}^{(r-1,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\!\!\!\!\!\!\!\!\!+\sum _{s\ge r} f_{0}^{(r-1,s)}({{\varvec{q}}})+\sum _{s\ge r} f_{1}^{(r-1,s)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})+\sum _{s\ge r} f_{2}^{(r-1,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}), \end{aligned} \end{aligned}$$

(19)

with $f_{l}^{(r-1, s)}\in \mathfrak {P}_{l,sK}\,$. By comparing formula (17) with (19), one immediately realizes that the assumption above is satisfied in the case with $r=1$ for what concerns the expansion of the initial Hamiltonian $\mathcal {H}^{(0)}$.

The r-th normalization step consists of three substeps, each of them involving a canonical transformation which is expressed in terms of the Lie series having $\chi _{0}^{(r)}$, $\chi _{1}^{(r)}$, $\chi _{2}^{(r)}$ as generating function, respectively. Therefore, the new Hamiltonian that is introduced at the end of the r-th normalization step is defined as follows:

$$\begin{aligned} \mathcal {H}^{(r)}=\exp \left( L_{\chi _2^{(r)}}\right) \exp \left( L_{\chi _1^{(r)}}\right) \exp \left( L_{\chi _0^{(r)}}\right) \mathcal {H}^{(r-1)}\, , \end{aligned}$$

(20)

where $\exp \left( L_{\chi } \right) \cdot = \sum _{j \ge 0 } (L_{\chi }^{j} \cdot )/j!$ is the Lie series operator, $L_{\chi } \cdot =\lbrace \cdot ,\chi \rbrace$ is the Lie derivative with respect to the dynamical function $\chi \,$, and $\lbrace \cdot ,\cdot \rbrace$ represents the Poisson bracket.

4.1.1 First substep (of the r-th normalization step)

The first substep aims to remove the term depending only on the angles¹²${{\varvec{q}}}$ up to trigonometric degree $rK\,$, i.e., $f_0^{(r-1,r)}$ (included in the first sum of the second row of (19)), which has to be considered as $\mathcal {O}(\varepsilon ^r)\,$. The first generating function $\chi _{0}^{(r)}({{\varvec{q}}})$ is determined by solving the following homological equation:

$$\begin{aligned} \lbrace {{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}},\chi _{0}^{(r)} \rbrace +f_{0}^{(r-1,r)}({{\varvec{q}}})=\left\langle f_{0}^{(r-1,r)} \right\rangle _{{{\varvec{q}}}}\,. \end{aligned}$$

(21)

Since $f_{0}^{(r-1,r)}\in \mathfrak {P}_{0,rK}\,$, its Fourier expansion can be written $f_{0}^{(r-1,r)}({{\varvec{q}}})=\sum _{|k|\le rK}c_{{{\varvec{k}}}}^{(r-1)} e^{i{{\varvec{k}}}\cdot {{\varvec{q}}}}$. Because of the homological Eq. (21), we find

$$\begin{aligned} \chi _{0}^{(r)}({{\varvec{q}}})=\sum _{0<|k|\le rK}\frac{c_{{{\varvec{k}}}}^{(r-1)}}{i\,{{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}} e^{i{{\varvec{k}}}\cdot {{\varvec{q}}}}\,; \end{aligned}$$

(22)

the above solution is well defined if the non-resonance condition

$$\begin{aligned} {{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}\ne 0 \quad \quad \forall \,0<{{\varvec{k}}}\le rK \end{aligned}$$

is satisfied. We can now apply the Lie series operator $\exp \big (L_{\chi _0^{(r)}}\big )$ to $\mathcal {H}^{(r-1)}\,$. This allows us to write the expansion of the new intermediate Hamiltonian as follows:

$$\begin{aligned} \begin{aligned}&\mathcal {H}^{(I;\,r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})=\exp \left( L_{\chi _0^{(r)}}\right) \mathcal {H}^{(r-1)}\\&\quad = \mathcal {E}^{(r)}+{{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(r-1)}\cdot {{\varvec{I}}}+\sum _{s\ge 0}\sum _{l\ge 3} f_{l}^{(I;\, r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\qquad +\sum _{s\ge r} f_{0}^{(I;\,r,s)}({{\varvec{q}}})+\sum _{s\ge r} f_{1}^{(I;\,r,s)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})+\sum _{s\ge r} f_{2}^{(I;\,r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,, \end{aligned} \end{aligned}$$

(23)

where (by abuse of notation) for the new canonical coordinates we adopt the same symbols as the old ones. From a practical point of view, the new Hamiltonian terms can be conveniently defined in such a way to mimic what is usually done in any programming language. First, we introduce the new summands as the old ones, so that $f_{l}^{(I;\,r,s)}=f_l^{(r-1,s)}$ $\forall \,l\ge 0\,$, $s\ge 0\,$. Hence, each term generated by Lie derivatives with respect to $\chi ^{(r)}_0$ is added to the corresponding class of functions. By a further abuse of notation, this is made by the following sequence¹³ of redefinitions:

$$\begin{aligned} f_{l-2j}^{(I;\,r,s+jr)} \hookleftarrow \frac{1}{j!} L_{\chi _{0}^{(r)}}^j f_{l}^{(r-1, s)} \qquad \forall \,l\ge 0 ,\,1\le j\le \lfloor l/2 \rfloor ,\,s\ge 0\, , \end{aligned}$$

(24)

where with the notation $a \hookleftarrow b$ we mean that the quantity a is redefined so as to be equal $a + b\,$. In fact, since $\chi _0^{(r)}\in \mathfrak {P}_{0,rK}\,$, Lemma 4.1 ensures that each application of the Lie derivative operator $L_{\chi _{0}^{(r)}}$ decreases by 1 the degree in ${{\varvec{p}}}$ (that is obviously equivalent to 2 in the square root of the actions), while the trigonometrical degree in the angles ${{\varvec{q}}}$ is increased by $rK\,$. By using repeatedly such a simple rule, one can easily verify that $f_{l}^{(I;\,r,s)}\in \mathfrak {P}_{l,sK}$ $\forall \, l\ge 0,\,s\ge 0\,$. Moreover, due to the homological Eq. (21), we set $f_{0}^{(I;\,r,r)}=0$ and update the energy value in such a way that $\mathcal {E}^{(r)}=\mathcal {E}^{(r-1)}+\left\langle f_{0}^{(r-1,r)} \right\rangle _{{{\varvec{q}}}}\,$, where $\left\langle \cdot \right\rangle _{{{\varvec{q}}}}$ is used to denote the angular average with respect to ${{\varvec{q}}}$.

4.1.2 Second substep (of the r-th normalization step)

The second substep aims to remove the term that is linear in $\sqrt{{{\varvec{I}}}}$ and independent on ${{\varvec{p}}}\,$, i.e., $f_1^{(I; r,r)}$, which is included in the second sum appearing in the second row of (23). The second generating function $\chi _{1}^{(r)}({{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}})$ is determined solving the following homological¹⁴ equation:

$$\begin{aligned} \lbrace {{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}}+ {{\varvec{\Omega }}}^{(r-1)}\cdot {{\varvec{I}}},\chi _{1}^{(r)} \rbrace +f_{1}^{(I;\,r,r)}({{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}})=0\,. \end{aligned}$$

(25)

Since $f_{1}^{(I;\,r,r)}\in \mathfrak {P}_{1,rK}\,$, we can write its expansion as

$$\begin{aligned} f_{1}^{(I;\,r,r)}({{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}})= \sum _{ 0\le |{{\varvec{k}}}|\le rK-1}\ \sum _{j=1}^{n_2} \sqrt{I_j}\left[ c_{{{\varvec{k}}},j}^{(+)}e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}+\alpha _j\right) } +c_{{{\varvec{k}}},j}^{(-)}e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}-\alpha _j\right) }\right] \,; \end{aligned}$$

due to the homological Eq. (25), we find

$$\begin{aligned} \chi _{1}^{(r)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}){} & {} = \sum _{ 0\le {{\varvec{k}}}\le rK-1}\sum _{j=1}^{n_2}\sqrt{I_j}\Bigg [\frac{c_{{{\varvec{k}}},j}^{(+)}}{i\,\left( {{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}+\Omega _j^{(r-1)}\right) }e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}+\alpha _j\right) } \nonumber \\{} & {} \quad +\frac{c_{{{\varvec{k}}},j}^{(-)}}{i\,\left( {{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}-\Omega _j^{(r-1)}\right) }e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}-\alpha _j\right) }\Bigg ]\,. \end{aligned}$$

(26)

The above expression is well defined provided that the first Melnikov non-resonance condition is satisfied, i.e.,

$$\begin{aligned} \min _{\begin{array}{c} 0<|{{\varvec{k}}}|\le rK-1 \\ |{{\varvec{l}}}|=1 \end{array}} \left| {{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}+{{\varvec{l}}}\cdot {{\varvec{\Omega }}}^{(r-1)}\right| \ge \frac{\gamma }{(rK)^\tau } \quad \quad \textrm{and}\quad \quad \min _{|{{\varvec{l}}}|=1}\left| \,{{\varvec{l}}}\cdot {{\varvec{\Omega }}}^{(r-1)}\right| \ge \gamma \,, \end{aligned}$$

(27)

for a pair of fixed values of $\gamma > 0$ and $\tau >n_1-1\,$ (see Locatelli et al. (2022) and reference therein).

We can now apply the transformation $\exp \left( L_{\chi _1^{(r)}}\right)$ to the Hamiltonian $\mathcal {H}^{(I;\,r)}\,$. By the usual abuse of notation (i.e., the new canonical coordinates are denoted with the same symbols of the old ones), the expansion of the new Hamiltonian can be written as

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(II;\,r)}&({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})=\exp \left( L_{\chi _1^{(r)}}\right) \mathcal {H}^{(I;\,r)}\\&\quad = \mathcal {E}^{(r)}+{{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(r-1)}\cdot {{\varvec{I}}}+\sum _{s\ge 0}\sum _{l\ge 3} f_{l}^{(II;\, r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\qquad +\sum _{s\ge r+1} f_{0}^{(II;\,r,s)}({{\varvec{q}}})+\sum _{s\ge r} f_{1}^{(II;\,r,s)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})+\sum _{s\ge r} f_{2}^{(II;\,\,r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,, \end{aligned} \end{aligned}$$

(28)

where in the last row of the previous formula, it is possible to start the first sum from $r+1$ instead of $r\,$, being $f_{0}^{(II;\,r,r)}=f_{0}^{(I;\,r,r)}=0\,$. In an analogous way as in the first substep, it is convenient to first define the new Hamiltonian terms as the old ones, i.e., $f_{l}^{(II;\,r,s)}=f_l^{(I;\,r,s)}$ $\forall \,l\ge 0\,$, $s\ge 0\,$. Hence, each term generated by the Lie derivatives with respect to $\chi ^{(r)}_1$ is added to the corresponding class of functions. This is made by the following sequence¹⁵ of redefinitions:

$$\begin{aligned} \begin{aligned}&f_{l-j}^{(II;\,r,s+jr)} \hookleftarrow \frac{1}{j!} L_{\chi _{1}^{(r)}}^j f_{l}^{(I;\,r, s)} \qquad \forall \,l\ge 0,\,1\le j\le l,\,s\ge 0\,,\\&f_{0}^{(II;\,r,2r)} \hookleftarrow \frac{1}{2} L_{\chi _{1}^{(r)}}^2 \left( {{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}}+\Omega ^{(r-1)}\cdot {{\varvec{I}}}\right) \,. \end{aligned} \end{aligned}$$

(29)

In fact, since $\chi _1^{(r)}\in \mathfrak {P}_{1,rK}$ is linear in $\sqrt{{{\varvec{I}}}}$, each application of the Lie derivative operator $L_{\chi _{1}^{(r)}}$ decreases by 1 the degree in square root of the actions, while the trigonometrical degree in the angles is increased by $rK\,$; such a rule holds because of Lemma 4.1. Moreover, thanks to the homological Eq. (25), one can easily remark that $f_{1}^{(II;\,r,r)}=0\,$. A repeated application of Lemma 4.1 allows us to verify also that $f_{l}^{(II;\,r,s)}\in \mathfrak {P}_{l,sK}\,$, $\forall l\ge 0,\,s\ge 0\,$.

4.1.3 Third substep (of the r-th normalization step)

The last substep aims to remove the term $f_{2}^{(II;\, r,r)}$ which is quadratic in the square root of the actions (i.e., either quadratic in $\sqrt{{{\varvec{I}}}}$ or linear in ${{\varvec{p}}}\,$) and included in the third sum appearing in the second row of (28). The third generating function $\chi _{2}^{(r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}})$ is determined by solving the following homological equation:

$$\begin{aligned} \lbrace {{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}}+ {{\varvec{\Omega }}}^{(r-1)}\cdot {{\varvec{I}}},\chi _{2}^{(r)} \rbrace +f_{2}^{(II;\,r,r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}})=\left\langle f_{2}^{(II;\,r,r)} \right\rangle _{{{\varvec{q}}},{{\varvec{\alpha }}}}\,. \end{aligned}$$

(30)

Since $f_{2}^{(II;\,r,r)}\in \mathfrak {P}_{2,rK}\,$, we can write it (according to definition 4.1 with $2\widehat{m}+\widehat{l}=2$ and $s=r\,$) as follows:

$$\begin{aligned} f_{2}^{(II;\,r,r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}})&=\sum _{\begin{array}{c} {{\varvec{m}}}\in \mathbb {N}^{n_1}\\ |{{\varvec{m}}}|=1 \end{array}}\sum _{\begin{array}{c} {{\varvec{k}}}\in \mathbb {Z}^{n_1}\\ |{{\varvec{k}}}|\le rK \end{array}}c_{{{\varvec{m}}},{{\varvec{k}}}}\,\,{{\varvec{p}}}^{{{\varvec{m}}}}e^{i\,{{\varvec{k}}}\cdot {{\varvec{q}}}} \\&\quad +\sum _{\begin{array}{c} {{\varvec{l}}}\in \mathbb {N}^{n_2}\\ |{{\varvec{l}}}|=2 \end{array}}\sum _{\begin{array}{c} {{\varvec{k}}}\in \mathbb {Z}^{n_1}\\ |{{\varvec{k}}}|+|\widehat{{{\varvec{l}}}}|\le rK \end{array}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\sum _{\begin{array}{c} \quad \qquad \widehat{l}_j=-l_j,-l_j+2,\ldots ,l_j\\ \,\,\,\,\,j=1,\ldots ,n_2 \end{array}}\!\!\!\!\!\!\!\tilde{c}_{{{\varvec{l}}},{{\varvec{k}}},\widehat{{{\varvec{l}}}}}\,\,\left( \sqrt{{{\varvec{I}}}}\right) ^{{{\varvec{l}}}}e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}+\widehat{{{\varvec{l}}}}\cdot {{\varvec{\alpha }}}\right) }\, . \end{aligned}$$

Due to the homological Eq. (30), we obtain

$$\begin{aligned} \begin{aligned}&\chi _{2}^{(r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}}) =\sum _{\begin{array}{c} {{\varvec{m}}}\in \mathbb {N}^{n_1}\\ |{{\varvec{m}}}|=1 \end{array}}\sum _{\begin{array}{c} {{\varvec{k}}}\in \mathbb {Z}^{n_1}\\ 0<|{{\varvec{k}}}|\le rK \end{array}}\frac{c_{{{\varvec{m}}},{{\varvec{k}}}}\,\,{{\varvec{p}}}^{{{\varvec{m}}}}e^{i\,{{\varvec{k}}}\cdot {{\varvec{q}}}}}{i\,{{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}}\\&\quad +\sum _{\begin{array}{c} {{\varvec{l}}}\in \mathbb {N}^{n_2}\\ |{{\varvec{l}}}|=2 \end{array}}\sum _{\begin{array}{c} {{\varvec{k}}}\in \mathbb {Z}^{n_1}\\ 0<|{{\varvec{k}}}|+|\widehat{{{\varvec{l}}}}|\le rK \end{array}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\sum _{\begin{array}{c} \quad \qquad \widehat{l}_j=-l_j,-l_j+2,\ldots ,l_j\\ \,\,\,\,\,j=1,\ldots ,n_2 \end{array}}\frac{\tilde{c}_{{{\varvec{l}}},{{\varvec{k}}},\widehat{{{\varvec{l}}}}}\,\,\left( \sqrt{{{\varvec{I}}}}\right) ^{{{\varvec{l}}}}e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}+\widehat{{{\varvec{l}}}}\cdot {{\varvec{\alpha }}}\right) }}{i\left( {{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}+\widehat{{{\varvec{l}}}}\cdot {{\varvec{\Omega }}}^{(r-1)}\right) }\,, \end{aligned} \end{aligned}$$

(31)

provided that both the non-resonance condition and the Melnikov one of second kind are satisfied, i.e.,

$$\begin{aligned} {{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}\ne 0 \quad \forall \,0<{{\varvec{k}}}\le rK \,, \ \, \min _{\begin{array}{c} 0<|{{\varvec{k}}}|\le rK-2\\ |{{\varvec{l}}}|=2 \end{array}}\left| {{\varvec{k}}}\cdot {{\varvec{\omega }}}^{(r-1)}+{{\varvec{l}}}\cdot {{\varvec{\Omega }}}^{(r-1)} \right| \ge \frac{\gamma }{(rK)^\tau } \,, \end{aligned}$$

(32)

with the same values of the constant parameters $\gamma > 0$ and $\tau >n_1-1$ appearing in (27).

We can now apply the transformation $\exp \left( L_{\chi _2^{(r)}}\right)$ to the Hamiltonian $\mathcal {H}^{(II;\,r)}\,$. By the usual abuse of notation (i.e., the new canonical coordinates are denoted with the same symbols as the old ones), the expansion¹⁶ of the new Hamiltonian can be written as

$$\begin{aligned} \begin{aligned}&\mathcal {H}^{(r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})=\exp \left( L_{\chi _2^{(r)}}\right) \mathcal {H}^{(II;\,r)}\\&\quad = \mathcal {E}^{(r)}+{{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(r-1)}\cdot {{\varvec{I}}}+\sum _{s\ge 0}\sum _{l\ge 3} f_{l}^{(r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\qquad +\sum _{s\ge r+1} f_{0}^{(r,s)}({{\varvec{q}}})+\sum _{s\ge r+1} f_{1}^{(r,s)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})+\sum _{s\ge r} f_{2}^{(r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,. \end{aligned} \end{aligned}$$

(33)

Once again, it is convenient to first define the new Hamiltonian terms as the old ones, i.e., $f_{l}^{(r,s)}=f_l^{(II;\,r,s)}$ $\forall \,l\ge 0\,$, $s\ge 0\,$. Hence, each term generated by the Lie derivatives with respect to $\chi ^{(r)}_2$ is added to the corresponding class of functions. This is made by the following sequence¹⁷ of redefinitions:

$$\begin{aligned} \begin{aligned}&f_{l}^{(\,r,s+jr)} \hookleftarrow \frac{1}{j!} L_{\chi _{1}^{(r)}}^j f_{l}^{(II;\,r, s)} \qquad \forall \,l\ge 0,\, j\ge 1,\,s\ge 0\,,\\&f_{2}^{(r,jr)} \hookleftarrow \frac{1}{j!} L_{\chi _{2}^{(r)}}^j \left( {{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(r-1)}\cdot {{\varvec{I}}}\right) \qquad \forall \,j\ge 1\,. \end{aligned} \end{aligned}$$

(34)

In fact, since $\chi _2^{(r)}\in \mathfrak {P}_{2,rK}$ is either quadratic in $\sqrt{{{\varvec{I}}}}$ or linear in ${{\varvec{p}}}$, each application of the Lie derivative operator $L_{\chi _{2}^{(r)}}$ does not modify the degree in the square root of the actions, while the trigonometric degree in the angles is increased by $rK\,$. By applying Lemma 4.1, one can verify also that $f_{l}^{(r,s)}\in \mathfrak {P}_{l,sK}\,$, $\forall l\ge 0,\,s\ge 0\,$.

Because of the homological Eq. (30), it immediately follows that the term that cannot be removed, that is $f_{2}^{(r,r)}=\left\langle f_{2}^{(II;\,r,r)} \right\rangle _{{{\varvec{q}}},{{\varvec{\alpha }}}}\in \mathfrak {P}_{2,0}\,$, is exactly of the same type with respect to the main term that is linear in the actions, i.e., ${{\varvec{\omega }}}^{(r-1)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(r-1)}\cdot {{\varvec{I}}}$. It looks then natural to update the angular velocity vectors so that

$$\begin{aligned} {{\varvec{\omega }}}^{(r)}={{\varvec{\omega }}}^{(r-1)}+\nabla _{{{\varvec{p}}}}\left( \left\langle f_{2}^{(II;\,r,r)} \right\rangle _{{{\varvec{q}}},{{\varvec{\alpha }}}}\right) , \quad \, {{\varvec{\Omega }}}^{(r)}={{\varvec{\Omega }}}^{(r-1)}+\nabla _{{{\varvec{I}}}}\left( \left\langle f_{2}^{(II;\,r,r)} \right\rangle _{{{\varvec{q}}},{{\varvec{\alpha }}}}\right) , \end{aligned}$$

(35)

where, as usual, the symbols $\nabla _{{{\varvec{p}}}}$ and $\nabla _{{{\varvec{I}}}}$ denote the gradient with respect to the action variables ${{\varvec{p}}}$ and ${{\varvec{I}}}$, respectively, and to set $f_{2}^{(r,r)}=0\,$. Therefore, the expansion of the Hamiltonian $\mathcal {H}^{(r)}$ can be rewritten as

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})&= \mathcal {E}^{(r)}+{{\varvec{\omega }}}^{(r)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(r)}\cdot {{\varvec{I}}}+\sum _{s\ge 0}\sum _{l\ge 3} f_{l}^{(r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\quad +\sum _{s\ge r+1}\sum _{l=0}^2 f_{l}^{(r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,, \end{aligned} \end{aligned}$$

(36)

where $f_{l}^{(r, s)}\in \mathfrak {P}_{l,sK}\,$ and $\mathcal {E}^{(r)}\in \mathfrak {P}_{0,0}\,$ is a constant.

It is now possible to iterate the algorithm, by performing the (next) $(r+1)$-th normalization step. The convergence of this normal form algorithm is proved in Caracciolo (2022) under suitable conditions.

In order to implement such a kind of normalization algorithm with the aid of a computer, we have to deal with Hamiltonians including a finite number of summands in their expansions in Taylor–Fourier series. To fix the ideas, let us suppose that we set a truncation rule in such a way as to neglect all the terms with a trigonometric degree greater than $\mathcal {N}_SK$, for a fixed positive integer value of the parameter $\mathcal {N}_S\,$. After iteratively performing $\mathcal {N}_S$ steps of the constructive algorithm, we end up with an approximation of the Hamiltonian which is in the normal form corresponding to an elliptic torus, i.e.,

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(\mathcal {N}_S)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})&= \mathcal {E}^{(\mathcal {N}_S)}+{{\varvec{\omega }}}^{(\mathcal {N}_S)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(\mathcal {N}_S)}\cdot {{\varvec{I}}}+\sum _{s = 0}^{\mathcal {N}_S}\sum _{l\ge 3} f_{l}^{(\mathcal {N}_S,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\,. \end{aligned} \end{aligned}$$

(37)

The Hamiltonian $\mathcal {H}^{(\mathcal {N}_S)}$ represents the natural starting point for the application of a second (Birkhoff-like) algorithm, which aims to produce a new normal form in such a way to remove the dependence on the angles ${{\varvec{q}}}$, as explained in the next subsection.

4.2 Construction of the resonant normal form in such a way to average with respect to the angles $\textbf{q}$

Consider a Hamiltonian¹⁸$\mathcal {H}_B^{(0)}$ of the form:

$$\begin{aligned} \begin{aligned} \mathcal {H}_B^{(0)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})&= \mathcal {E}_B+{{\varvec{\omega }}}_B\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}_B\cdot {{\varvec{I}}}+\sum _{s= 0}^{\mathcal {N}_S}\sum _{l\ge 3} g_{l}^{(0,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\,, \end{aligned} \end{aligned}$$

(38)

where $\mathcal {E}_B$ is a constant term, representing the energy, $({{\varvec{p}}},{{\varvec{q}}})\in \mathbb {R}^{n_1}\times \mathbb {T}^{n_1}\,$, $({{\varvec{I}}},{{\varvec{\alpha }}})\in \mathbb {R}^{n_2}_{\ge 0}\times \mathbb {T}^{n_2}$ are action-angle variables, $({{\varvec{\omega }}}_B,{{\varvec{\Omega }}}_B)\in \mathbb {R}^{n_1}\times \mathbb {R}^{n_2}$ are the frequencies, $\mathcal {N}_S$ is a fixed positive integer (ruling the truncations in the Fourier series) and $g_{l}^{(0, s)}\in \mathfrak {P}_{l,sK}\,$, $\forall \,l\ge 0,\,0\le s\le \mathcal {N}_S\,$. In practice, we are starting from the normalized Hamiltonian of the previous subsection $\mathcal {H}^{(\mathcal {N}_S)}\,$, given by Eq. (37), where we have defined $\mathcal {H}^{(0)}_B:=\mathcal {H}^{(\mathcal {N}_S)}\,$, $\mathcal {E}_B:=\mathcal {E}^{(\mathcal {N}_S)}\,$, $({{\varvec{\omega }}}_B,{{\varvec{\Omega }}}_B):=({{\varvec{\omega }}}^{(\mathcal {N}_S)}, {{\varvec{\Omega }}}^{(\mathcal {N}_S)})$ and $g_{l}^{(0, s)}:=f_{l}^{(\mathcal {N}_S, s)}\in \mathfrak {P}_{l,sK}\,$ $\forall \,l\ge 0,\,0\le s\le \mathcal {N}_S\,$; this is done also in order to simplify the notation. By comparison with Eq. (37), it is easy to remark that $g_{l}^{(0, s)}:=f_{l}^{(\mathcal {N}_S, s)}=0\,$, $\forall \, 0\le l \le 2\,$, $1\le s \le \mathcal {N}_S\,$.

The aim of the present algorithm is to delete the dependence of $\mathcal {H}^{(0)}_B$ on the angles ${{\varvec{q}}}\,$, reducing by $n_1$ the number of degrees of freedom. In order to do this, we have to act on the terms $g_{l}^{(0,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})$ such that $s\ge 1$ and $l\ge 3$, removing their dependence on ${{\varvec{q}}}\,$; indeed, for $s=0\,$, the sum $\sum _{l\ge 3} g_{l}^{(0,0)}({{\varvec{p}}},{{\varvec{I}}})\,$ does not depend on the angles, thus it is already in normal form. This elimination can be done by a sequence of canonical transformations. If convergent, this would lead the Hamiltonian to the following final normal form:

$$\begin{aligned} \begin{aligned} \mathcal {H}_B^{(\infty )}(\tilde{{{\varvec{p}}}}, \tilde{{{\varvec{I}}}},\tilde{{{\varvec{\alpha }}}})&= \mathcal {E}_B+{{\varvec{\omega }}}_B\cdot \tilde{{{\varvec{p}}}}+{{\varvec{\Omega }}}_B\cdot \tilde{{{\varvec{I}}}}+ \sum _{s= 0}^{\infty }\sum _{l\ge 3} g_{l}^{(\infty ,s)}(\tilde{{{\varvec{p}}}},\tilde{{{\varvec{I}}}},\tilde{{{\varvec{\alpha }}}})\,, \end{aligned} \end{aligned}$$

(39)

where $(\tilde{{{\varvec{p}}}},\tilde{{{\varvec{I}}}},\tilde{{{\varvec{\alpha }}}})$ denote the new variables; it is evident that, having removed the dependence on $\tilde{{{\varvec{q}}}}\,$, the conjugate momenta vector $\tilde{{{\varvec{p}}}}$ is constant. However, as typical of the computational procedures à la Birkhoff, the constructive algorithm produces divergent series if the normalization is iterated infinitely many times. For this reason, it is convenient to look for an optimal order of normalization to which the algorithm is stopped. In our approach, we have not to consider such a problem, because we are dealing with truncated series; this is done in order to keep our discussion as close as possible to the practical implementations where the maximal degree in actions of the expansions is usually rather low.

The generic r-th step of this new normalization algorithm is defined as follows. After $r-1$ steps, the Hamiltonian (38) takes the form

$$\begin{aligned} \begin{aligned}&\mathcal {H}_B^{(r-1)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})= \mathcal {E}_B+{{\varvec{\omega }}}_B\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}_B\cdot {{\varvec{I}}}+\sum _{l\ge 3} g_{l}^{(r-1,0)}({{\varvec{p}}},{{\varvec{I}}})\\&\quad +\sum _{s=1}^{\mathcal {N}_S}\sum _{3\le l\le r+1} g_{l}^{(r-1,s)}({{\varvec{p}}},{{\varvec{I}}},{{\varvec{\alpha }}})+\sum _{s= 1}^{\mathcal {N}_S}\sum _{l\ge r+2} g_{l}^{(r-1,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,, \end{aligned} \end{aligned}$$

(40)

with $g_{l}^{(r-1,s)}\in \mathfrak {P}_{l,sK}\,$.

The r-th normalization step consists of a sequence of $\mathcal {N}_S$ substeps, each of them involving a canonical transformation which is expressed in terms of the Lie series having $\chi _B^{(j;\,r)}$ as generating function, with $j=1,\,\ldots ,\mathcal {N}_S\,$. Therefore, the new Hamiltonian introduced at the end of the r-th normalization step of this algorithm is defined as follows:

$$\begin{aligned} {\mathcal {H}}_{B}^{(r)}= \exp \left( L_{\chi _B^{(\mathcal {N}_S;\,r)}}\right) \ldots \exp \left( L_{\chi _B^{(3;\,r)}}\right) \exp \left( L_{\chi _B^{(2;\,r)}}\right) \exp \left( L_{\chi _B^{(1;\,r)}}\right) \mathcal {H}_B^{(r-1)}\,. \end{aligned}$$

(41)

The generating functions $\chi _{B}^{(j;\,r)}$ are defined so as to remove the dependence on ${{\varvec{q}}}$ from the perturbing term of lowest order in the square root of the actions, i.e., $\sum _{s= 1}^{\mathcal {N}_S} g_{r+2}^{(r-1,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\,$.

4.2.1 j-th substep of the r-th step of the algorithm constructing the resonant normal form

After $j-1$ substeps, the Hamiltonian can be written as follows:

$$\begin{aligned} \begin{aligned}&\mathcal {H}_B^{(j-1; \,r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})= \mathcal {E}_B+{{\varvec{\omega }}}_B\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}_B\cdot {{\varvec{I}}} +\sum _{l\ge 3} g_{l}^{(j-1;\,r,0)}({{\varvec{p}}},{{\varvec{I}}})\\&\quad +\sum _{s= 1}^{\mathcal {N}_S}\sum _{3\le l\le r+1} g_{l}^{(j-1;\,r,s)}({{\varvec{p}}},{{\varvec{I}}},{{\varvec{\alpha }}}) +\sum _{s=1}^{j-1} g_{r+2}^{(j-1;\,r,s)}({{\varvec{p}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\quad +\sum _{s= j}^{\mathcal {N}_S} g_{r+2}^{(j-1;\,r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})+ \sum _{s= 1}^{\mathcal {N}_S}\sum _{l\ge r+3} g_{l}^{(j-1;\,r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}), \end{aligned} \end{aligned}$$

(42)

where, for $j=1\,$, we set $\mathcal {H}_{B}^{(0;\,r)}:=\mathcal {H}_{B}^{(r-1)}\,$ and $g_{l}^{(0;\,r,s)}=g_{l}^{(r-1,s)}\,$, $\forall l \ge 0$, $\forall \, 0\le s\le \mathcal {N}_S\,$.

The j-th substep generating function $\chi _{B}^{(j;\,r)}$ is determined by the following homological equation:

$$\begin{aligned} \begin{aligned} \lbrace {{\varvec{\omega }}}_{B}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}_{B}\cdot {{\varvec{I}}}\,,\chi _{B}^{(j;\,r)} \rbrace +g_{r+2}^{(j-1;\,r,j)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})= \left\langle g_{r+2}^{(j-1;\,r,j)} \right\rangle _{{{\varvec{q}}}} \end{aligned}\,. \end{aligned}$$

(43)

Proceeding in a similar way as in the description of the third substep of Sect. 4.1, first we write the expansion of the perturbing function in the form

$$\begin{aligned}{} & {} g_{r+2}^{(j-1; \,r,j)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}}) = \sum _{2|{{\varvec{m}}}|+|{{\varvec{l}}}|=r+2}\,\, \sum _{{{\varvec{m}}}\in \mathbb {N}^{n_1}}\, \sum _{{{\varvec{l}}}\in \mathbb {N}^{n_2}}\, \sum _{\begin{array}{c} {{\varvec{k}}}\in \mathbb {Z}^{n_1}\\ |{{\varvec{k}}}|+|\widehat{{{\varvec{l}}}}|\le jK \end{array}} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \nonumber \\{} & {} \quad \sum _{\begin{array}{c} \quad \qquad \widehat{l}_{j_2}=-l_{j_2},-l_{j_2}+2,\ldots ,l_{j_2}\\ \,\,\,\,\,j_2=1,\ldots ,n_2 \end{array}} \!\!\!\!\!\!\!c_{{{\varvec{m}}},{{\varvec{l}}},{{\varvec{k}}},\widehat{{{\varvec{l}}}}}\,\, {{\varvec{p}}}^{{{\varvec{m}}}}\left( \sqrt{{{\varvec{I}}}}\right) ^{{{\varvec{l}}}} e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}+\widehat{{{\varvec{l}}}}\cdot {{\varvec{\alpha }}}\right) }\,. \end{aligned}$$

(44)

The solution of the homological Eq. (43) is then

$$\begin{aligned} \begin{aligned}&\chi _{B}^{(j;\, r)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}}, {{\varvec{\alpha }}}) =\sum _{2|{{\varvec{m}}}|+|{{\varvec{l}}}|=r+2}\,\, \sum _{{{\varvec{m}}}\in \mathbb {N}^{n_1}}\, \sum _{{{\varvec{l}}}\in \mathbb {N}^{n_2}}\, \sum _{\begin{array}{c} {{\varvec{k}}}\in \mathbb {Z}^{n_1},\,|{{\varvec{k}}}|>0\\ |{{\varvec{k}}}|+|\widehat{{{\varvec{l}}}}|\le jK \end{array}} \\&\ \sum _{\begin{array}{c} \quad \qquad \widehat{l}_{j_2}=-l_{j_2},-l_{j_2}+2,\ldots ,l_{j_2}\\ \,\,\,j_2=1,\ldots ,n_2 \end{array}} \!\!\!\!\!\!\!\,\,\,\,\, \frac{c_{{{\varvec{m}}},{{\varvec{l}}},{{\varvec{k}}},\widehat{{{\varvec{l}}}}}}{i\left( {{\varvec{k}}}\cdot {{\varvec{\omega }}}_B+\widehat{{{\varvec{l}}}}\cdot {{\varvec{\Omega }}}_B \right) } \,\,\,\, {{\varvec{p}}}^{{{\varvec{m}}}}\left( \sqrt{{{\varvec{I}}}}\right) ^{{{\varvec{l}}}} e^{i\left( {{\varvec{k}}}\cdot {{\varvec{q}}}+\widehat{{{\varvec{l}}}}\cdot {{\varvec{\alpha }}}\right) }\,. \end{aligned} \end{aligned}$$

(45)

We can now apply the transformation $\exp \big (L_{\chi _B^{(j;\,r)}}\big )$ to the Hamiltonian $\mathcal {H}_B^{(j-1;\,r)}\,$. By the usual abuse of notation (i.e., the new canonical coordinates are denoted with the same symbols of the old ones), the expansion of the new Hamiltonian can be written as

$$\begin{aligned} \begin{aligned} \mathcal {H}_B^{(j;\,r)}&({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})= \exp \left( L_{\chi _B^{(j;\,r)}}\right) \mathcal {H}_B^{(j-1;\,r)}\\&\quad =\mathcal {E}_B+{{\varvec{\omega }}}_B\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}_B\cdot {{\varvec{I}}} +\sum _{l\ge 3} g_{l}^{(j;\,r,0)}({{\varvec{p}}},{{\varvec{I}}})\\&\qquad +\sum _{s=1}^{\mathcal {N}_S}\sum _{3\le l\le r+1} g_{l}^{(j;\,r,s)}({{\varvec{p}}},{{\varvec{I}}},{{\varvec{\alpha }}}) +\sum _{s=1}^{j} g_{r+2}^{(j;\,r,s)}({{\varvec{p}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\qquad +\sum _{s= j+1}^{\mathcal {N}_S} g_{r+2}^{(j;\,r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) +\sum _{s= 1}^{\mathcal {N}_S}\sum _{l\ge r+3} g_{l}^{(j;\,r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,. \end{aligned} \end{aligned}$$

(46)

In a similar way to what has been done previously, it is convenient to first define the new Hamiltonian terms as the old ones, i.e., $g_{l}^{(j;\,r,s)}=g_l^{(j-1;\,r,s)}$ $\forall \,l\ge 0\,$, $0\le s\le \mathcal {N}_S\,$; hence, each term generated by the Lie derivatives with respect to $\chi _B^{(j;\,r)}$ is added to the corresponding class of functions. This is made by the following sequence¹⁹ of redefinitions

$$\begin{aligned} \begin{aligned}&g_{l+mr}^{(j;\,r,s+mj)} \hookleftarrow \frac{1}{m!} L_{\chi _{B}^{(j;\,r)}}^m g_{l}^{(j;\,r, s)} \qquad \forall \,l\ge 0,\, 1\le m\le \lfloor (\mathcal {N}_S-s)/j \rfloor , \,0\le s\le \mathcal {N}_S\,,\\&\quad g_{2+mr}^{(j;\,r,mj)} \hookleftarrow \frac{1}{m!} L_{\chi _{B}^{(j;\,r)}}^m \left( {{\varvec{\omega }}}_B\cdot {{\varvec{p}}}+\Omega _B\cdot {{\varvec{I}}}\right) \qquad \forall \, 1\le m\le \lfloor \mathcal {N}_S/j\rfloor \,. \end{aligned} \end{aligned}$$

(47)

In fact, since $\chi _B^{(j;\,r)}\in \mathfrak {P}_{r+2,jK}\,$, each application of the Lie derivative operator $L_{\chi _B^{(j;\,r)}}$ increases the degree in square root of the actions and the trigonometrical degree in the angles by r and $jK\,$, respectively. Moreover, thanks to the homological Eq. (43) and the second rule included in formula (47) (in the case with $m=1$), one can easily remark that $g_{r+2}^{(j;\,r,j)}=\left\langle g_{r+2}^{(j-1;\,r,j)} \right\rangle _{{{\varvec{q}}}}\,$. By applying Lemma 4.1 one can verify also that $g_{l}^{(j;\,r,s)}\in \mathfrak {P}_{l,sK}\,$, $\forall l\ge 0,\,s\ge 0\,$.

The r-th step of the algorithm constructing the resonant normal form is completed at the end of the iterative repetition of the j-th substep for $j=1,\,\ldots ,\,\mathcal {N}_S\,$. Therefore, the expansion of the Hamiltonian can be written in the following form:

$$\begin{aligned} \begin{aligned} {\mathcal {H}}_B^{(r)}&({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})= \exp \left( L_{\chi _B^{(\mathcal {N}_S;\,r)}}\right) \cdots \exp \left( L_{\chi _B^{(1;\,r)}}\right) \mathcal {H}_B^{(r-1)}\\&\quad = \mathcal {E}_B+{{\varvec{\omega }}}_B\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}_B\cdot {{\varvec{I}}} +\sum _{l\ge 3} g_{l}^{(r,0)}({{\varvec{p}}},{{\varvec{I}}})\\&\qquad +\sum _{s= 1}^{\mathcal {N}_S}\sum _{3\le l\le r+2} g_{l}^{(r,s)}({{\varvec{p}}},{{\varvec{I}}},{{\varvec{\alpha }}}) +\sum _{s= 1}^{\mathcal {N}_S}\sum _{l\ge r+3} g_{l}^{(r,s)}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,, \end{aligned} \end{aligned}$$

(48)

where $g_{l}^{(r,s)}:=g_{l}^{(\mathcal {N}_S;\,r,s)}\,$, $\forall \,l\ge 0\,$, $0\le s\le \mathcal {N}_S\,$. Then, the normalization algorithm can be iteratively repeated. Since we are interested in the computer implementation, we consider finite sequences of Hamiltonians whose expansion is truncated up to a finite degree, say, $\mathcal {N}_L$ in the square root of the actions. Therefore, the iteration of $\mathcal {N}_L-2$ normalization steps of the algorithm constructing the resonant normal form are sufficient to obtain

$$\begin{aligned} \mathcal {H}_{B}^{(\mathcal {N}_L-2)}({{\varvec{p}}}, {{\varvec{I}}},{{\varvec{\alpha }}}) = \mathcal {E}_B+{{\varvec{\omega }}}_B\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}_B\cdot {{\varvec{I}}}+ \sum _{s=0}^{\mathcal {N}_S}\sum _{l=3}^{\mathcal {N}_L} g_{l}^{(\mathcal {N}_L-2,s)}({{\varvec{p}}},{{\varvec{I}}},{{\varvec{\alpha }}})\,. \end{aligned}$$

(49)

The Hamiltonian (49) does not depend on the angles ${{\varvec{q}}}$. Therefore, the corresponding actions ${{\varvec{p}}}$ are constant and they can be considered as parameters whose values are fixed by the initial conditions; this allows us to decrease the number of degrees of freedom by $n_1\,$, passing from $n_1+n_2$ to $n_2$.

5 Application of the normalization algorithms to the secular quasi-periodic restricted model of the dynamics of $\upsilon$-And $\textbf{b}$

The SQPR model can be reformulated in such a way as to resume the form of a Hamiltonian of the type (17), to which we can sequentially apply both normalization procedures described in the two previous subsections. In fact, the canonical change of variables

$$\begin{aligned} \begin{aligned}&\xi _1=\sqrt{2 I_1}\cos (\alpha _1) \,, \qquad{} & {} \eta _1=\sqrt{2 I_1}\sin (\alpha _1)\,,\\&P_1=\sqrt{2 I_2}\cos (\alpha _2)\,, \qquad{} & {} Q_1=\sqrt{2 I_2}\sin (\alpha _2)\,, \end{aligned} \end{aligned}$$

(50)

allows to rewrite the expansion of the SQPR Hamiltonian (13) as follows:

$$\begin{aligned} \begin{aligned}&\mathcal {H}_{sec,\, 2+3/2}({{\varvec{p}}},{{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) = \omega _3\,p_3 + \omega _4\,p_4 + \omega _5\,p_5 \\&+\!\!\!\!\!\!\!\!\ \sum _{\begin{array}{c} l_1+l_2=0\\ (l_1,\,l_2)\in \mathbb {N}^2 \end{array}}^{\mathcal {N}_L} \sum _{\begin{array}{c} (k_3,\,k_4,\,k_5)\in \mathbb {Z}^3\\ |{{\varvec{k}}}|\le \mathcal {N}_SK \end{array}} \!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\sum _{\begin{array}{c} \qquad \qquad k_j=-l_j,-l_j+2,\ldots ,l_j\\ j=1,\,2 \end{array}} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! c_{{{\varvec{l}}},{{\varvec{k}}}} (\sqrt{I_1})^{l_1}(\sqrt{I_2})^{l_2} e^{i(k_1\alpha _1+k_2\alpha _2+k_3 q_3 + k_4 q_4 + k_5 q_5)}, \end{aligned} \end{aligned}$$

(51)

where ${{\varvec{k}}}=(k_1,\ldots ,k_5)\in \mathbb {Z}^5\,$. The r.h.s. of the above equation can be expressed in the general and more compact form described in Eq. (17), by setting $n_1=3\,$, $n_2=2\,$, ${{\varvec{\omega }}}^{(0)}={{\varvec{\omega }}}=(\omega _3,\omega _4,\omega _5)\in \mathbb {R}^3$, that are the fundamental frequencies of the two outer planets (described in Eq. (8)), while ${{\varvec{\Omega }}}^{(0)}\in \mathbb {R}^2$ can be easily determined by performing the so-called diagonalization of the Hamiltonian part that is quadratic in the square root of the actions ${{\varvec{I}}}$ and not depending on the angles ${{\varvec{q}}}$ (see, e.g., Giorgilli et al. 1989). In the equation above, the parameters $\mathcal {N}_L$ and $\mathcal {N}_S$ define the truncation order of the expansions in Taylor and Fourier series, respectively, in such a way to represent on the computer just a finite number of terms that are not too many to handle with; in our computations we fix $\mathcal {N}_L=6$ as maximal power degree in square root of the actions and we include Fourier terms up to a maximal trigonometric degree of 8, putting $\mathcal {N}_S=4\,$, $K=2\,$. We recall that setting $K=2$ is quite natural for Hamiltonian systems close to stable equilibria as it is for models describing the secular planetary dynamics, see, e.g., Giorgilli et al. (2017). Let us also remark that a simple reordering of the summands according to the total trigonometric degree $|{{\varvec{k}}}|$ in the angles $({{\varvec{q}}},{{\varvec{\alpha }}})$ allows us to represent the second row of formula (51) as a sum of Hamiltonian terms each of them is belonging to a functions class of type $\mathfrak {P}_{l,sK}\,$, which is unique for any positive integer K if we ask for the minimality of the index s. These comments can be used all together in order to formally verify that the new expansion of $\mathcal {H}_{sec,\, 2+3/2}$ in (51) can be finally reexpressed in the same form as $\mathcal {H}^{(0)}$ in (17).

Furthermore, in the case of our SQPR model of the secular dynamics of $\upsilon$-And $\textbf{b}$, the only term depending on the action variables ${{\varvec{p}}}$ (that are the so-called dummy variables) is ${{\varvec{\omega }}}^{(0)}\cdot {{\varvec{p}}}\,$; thus, none of the Hamiltonian term $f_{l}^{(0,s)}$ depends on ${{\varvec{p}}}\,$. This fact would allow to introduce some simplification in the computational algorithm. For instance, the value of the angular velocity vector ${{\varvec{\omega }}}^{(0)}$ is not modified during the first normalization procedure (i.e., the algorithmic construction of the elliptic tori) and it remains equal to its initial value, given by the fundamental frequencies described in (8). Therefore, the expansion of the starting Hamiltonian in the special case of our SQPR model can be rewritten as

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(0)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})&= \mathcal {E}^{(0)}+{{\varvec{\omega }}}^{(0)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(0)}\cdot {{\varvec{I}}}+\sum _{s= 0}^{\mathcal {N}_S}\sum _{l= 3}^{\mathcal {N}_L} f_{l}^{(0,s)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\\&\quad +\sum _{s= 1}^{\mathcal {N}_S}\sum _{l=0}^{2} f_{l}^{(0,s)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}}) \,; \end{aligned} \end{aligned}$$

(52)

however, in our opinion, for what concerns the general description of the previous subsections it has been worth to consider also an eventual dependence of $f_{l}^{(0,s)}$ on ${{\varvec{p}}}\,$ in order to keep the discussion of the constructive procedure as general as possible.

The first algorithm to be applied aims to construct the normal form corresponding to an invariant elliptic torus. It starts from the Hamiltonian $\mathcal {H}_{sec,\, 2+3/2}$ rewritten in the same form as $\mathcal {H}^{(0)}$ in (17) (more precisely as in (52)) and its computational procedure is fully detailed in Sect. 4.1. Therefore, we perform $\mathcal {N}_S$ normalization steps of this first normalization algorithm. This allows us to bring the Hamiltonian in the following (intermediate) normal form:

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(\mathcal {N}_S)}({{\varvec{p}}},{{\varvec{q}}}, {{\varvec{I}}},{{\varvec{\alpha }}})&= \mathcal {E}^{(\mathcal {N}_S)}+{{\varvec{\omega }}}^{(\mathcal {N}_S)}\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}^{(\mathcal {N}_S)}\cdot {{\varvec{I}}}+\sum _{s= 0}^{\mathcal {N}_S}\sum _{l= 3}^{\mathcal {N}_L} f_{l}^{(\mathcal {N}_S,s)}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\, , \end{aligned} \end{aligned}$$

where $f_{l}^{(\mathcal {N}_S,s)}\in \mathfrak {P}_{l,sK}$ $\forall \ l=3,\,\ldots \,,\,\mathcal {N}_L,\>s=0\,,\,\ldots \,,\,\mathcal {N}_S$ and the angular velocity vector related to the angles ${{\varvec{q}}}$ is such that ${{\varvec{\omega }}}^{(\mathcal {N}_S)}={{\varvec{\omega }}}^{(0)}=(\omega _3, \omega _4, \omega _5)$, whose components are given in (8).

It is now possible to apply the second algorithm aiming to construct a resonant normal form where the dependence on the angles ${{\varvec{q}}}=(q_3, q_4, q_5)$ is completely removed. Such a normalization starts from the Hamiltonian $\mathcal {H}^{(\mathcal {N}_S)}$ obtained after the first normalization procedure. Therefore, we perform $\mathcal {N}_L-2$ normalization steps of the above algorithm, each of them involving $\mathcal {N}_S$ substeps as described in Sect. 4.2; this allows us to bring the Hamiltonian in the following (final) normal form:

$$\begin{aligned} \begin{aligned} \mathcal {H}_{2DOF}({{\varvec{p}}}, {{\varvec{I}}},{{\varvec{\alpha }}})&= \mathcal {E}_B+{{\varvec{\omega }}}_B\cdot {{\varvec{p}}}+{{\varvec{\Omega }}}_B\cdot {{\varvec{I}}}\\&\quad +\sum _{l= 3}^{\mathcal {N}_L} g_{l}^{(\mathcal {N}_L-2,0)}({{\varvec{I}}})+\sum _{s= 1}^{\mathcal {N}_S}\sum _{l=3}^{\mathcal {N}_L} g_{l}^{(\mathcal {N}_L-2,s)}({{\varvec{I}}},{{\varvec{\alpha }}})\,, \end{aligned} \end{aligned}$$

(53)

where $g_{l}^{(\mathcal {N}_L-2,s)}\in \mathfrak {P}_{l,sK}$ $\forall \ l=3,\,\ldots \,,\,\mathcal {N}_L,\>s=0\,,\,\ldots \,,\,\mathcal {N}_S$ and it still holds true that ${{\varvec{\omega }}}_B={{\varvec{\omega }}}^{(0)}\,$.

All the algebraic manipulations that are prescribed by the normal form algorithms have been performed by using the symbolic manipulator Mathematica as a programming framework.

We emphasize that $\mathcal {H}_{2DOF}$ is an integrable Hamiltonian. In fact, due to the preservation of the total angular momentum, discussed in Remark 3.1, the following invariance law²⁰ is satisfied:

$$\begin{aligned} \frac{\partial \,\mathcal {H}_{2DOF}}{\partial \alpha _1}+ \frac{\partial \,\mathcal {H}_{2DOF}}{\partial \alpha _2}=0; \end{aligned}$$

(54)

thus, from the Hamilton’s equations for $\mathcal {H}_{2DOF}\,$, we can immediately deduce that $I_1+I_2$ is a constant of motion. Therefore, $\mathcal {H}_{2DOF}$ is integrable because of the Liouville theorem (see, e.g., Giorgilli 2022), since it admits a complete system of constants of motion in involution, that are the dummy variables ${{\varvec{p}}}$ (which could be disregarded in (53), reducing the model to 2 DOF), $I_1+I_2$ and the Hamiltonian itself.

In view of the numerical explorations of the dynamical evolution of our new model described by the integrable Hamiltonian $\mathcal {H}_{2DOF}({{\varvec{p}}},{{\varvec{I}}},{{\varvec{\alpha }}})$, it is convenient to introduce the canonical transformations related to the so-called semi-analytic method of integration for the equations of motion (see, e.g., Giorgilli et al. 2017). In order to fix the ideas, let us focus on the second algorithm, designed to construct a resonant normal form. This normalization procedure can be summarized by the transformation that is obtained by iteratively applying all the Lie series to the canonical variables. This is done as follows:

$$\begin{aligned} \begin{aligned} I_i&= \exp \left( L_{\chi _B^{(\mathcal {N}_S;\,\mathcal {N}_L-2)}}\right) \dots \exp \left( L_{\chi _B^{(2;\,\mathcal {N}_L-2)}}\right) \exp \left( L_{\chi _B^{(1;\,\mathcal {N}_L-2)}}\right) \ldots \dots \\&\quad \ldots \exp \left( L_{\chi _B^{(\mathcal {N}_S;\,1)}}\right) \dots \exp \left( L_{\chi _B^{(2;\,1)}}\right) \exp \left( L_{\chi _B^{(1;\,1)}}\right) I_i\bigg |_{\begin{array}{c} {{{\varvec{I}}}=\tilde{{{\varvec{I}}}}} \\ {{{\varvec{\alpha }}}=\tilde{{{\varvec{\alpha }}}}} \end{array}}\,,\\ \alpha _i&= \exp \left( L_{\chi _B^{(\mathcal {N}_S;\,\mathcal {N}_L-2)}}\right) \dots \exp \left( L_{\chi _B^{(2;\,\mathcal {N}_L-2)}}\right) \exp \left( L_{\chi _B^{(1;\,\mathcal {N}_L-2)}}\right) \ldots \dots \\&\quad \ldots \exp \left( L_{\chi _B^{(\mathcal {N}_S;\,1)}}\right) \dots \exp \left( L_{\chi _B^{(2;\,1)}}\right) \exp \left( L_{\chi _B^{(1;\,1)}}\right) \alpha _i\bigg |_{\begin{array}{c} {{{\varvec{I}}}=\tilde{{{\varvec{I}}}}} \\ {{{\varvec{\alpha }}}=\tilde{{{\varvec{\alpha }}}}} \end{array}}\,, \end{aligned} \end{aligned}$$

(55)

for $i=1,2\,$. We introduce the symbol $\mathscr {C}_B$ to denote the change of coordinates²¹ defined by the above expressions, i.e., $({{\varvec{I}}},{{\varvec{\alpha }}})=\mathscr {C}_B({{\varvec{q}}},\tilde{{{\varvec{I}}}},\tilde{{{\varvec{\alpha }}}})\,$. We can proceed in the same way for what concerns the algorithm constructing the normal form corresponding to an invariant elliptic torus. In fact, we first introduce the application of all the Lie series to the canonical variables in such a way to write, $\forall \ i=1,2\,$,

$$\begin{aligned} \begin{aligned} I_i=&\exp \left( L_{\chi _2^{(\mathcal {N}_S)}}\right) \exp \left( L_{\chi _1^{(\mathcal {N}_S)}}\right) \exp \left( L_{\chi _0^{(\mathcal {N}_S)}}\right) \ldots \\&\quad \ldots \exp \left( L_{\chi _2^{(1)}}\right) \exp \left( L_{\chi _1^{(1)}}\right) \exp \left( L_{\chi _0^{(1)}}\right) I_i\bigg |_{\begin{array}{c} {{{\varvec{I}}}=\hat{{{\varvec{I}}}}} \\ {{{\varvec{\alpha }}}=\hat{{{\varvec{\alpha }}}}} \end{array}}\,,\\ \alpha _i=&\exp \left( L_{\chi _2^{(\mathcal {N}_S)}}\right) \exp \left( L_{\chi _1^{(\mathcal {N}_S)}}\right) \exp \left( L_{\chi _0^{(\mathcal {N}_S)}}\right) \ldots \\&\quad \ldots \exp \left( L_{\chi _2^{(1)}}\right) \exp \left( L_{\chi _1^{(1)}}\right) \exp \left( L_{\chi _0^{(1)}}\right) \alpha _i\bigg |_{\begin{array}{c} {{{\varvec{I}}}=\hat{{{\varvec{I}}}}} \\ {{{\varvec{\alpha }}}=\hat{{{\varvec{\alpha }}}}} \end{array}}\,; \end{aligned} \end{aligned}$$

(56)

finally, we use the symbol $\mathscr {C}$ to summarize the whole change of coordinates that is defined by the whole expression above, i.e., $({{\varvec{I}}},{{\varvec{\alpha }}})=\mathscr {C}({{\varvec{q}}},\hat{{{\varvec{I}}}},\hat{{{\varvec{\alpha }}}})\,$. Let us now introduce the function $\mathscr {F}:\mathbb {T}^3\times \mathbb {R}_{\ge 0}^{2}\times \mathbb {T}^{2} \rightarrow \mathbb {R}_{\ge 0}^{2}\times \mathbb {T}^{2}$, which is defined so that

$$\begin{aligned} \mathscr {F}({{\varvec{q}}},\tilde{{{\varvec{I}}}},\tilde{{{\varvec{\alpha }}}})= \mathscr {C}\big ({{\varvec{q}}},\mathscr {C}_B({{\varvec{q}}},\tilde{{{\varvec{I}}}},\tilde{{{\varvec{\alpha }}}})\big ), \end{aligned}$$

(57)

where we have omitted to put the $\tilde{}$ symbol on top of ${{\varvec{q}}}$ in order to stress that the angles ${{\varvec{q}}}$ are not affected by the change of coordinates. Moreover, let also introduce the symbol ${{\mathcal {A}}}$ to denote the usual canonical transformation defining the action-angle variables for the harmonic oscillator, i.e., by formula (50), in our case this means that

$$\begin{aligned} {{\mathcal {A}}}({{\varvec{I}}},{{\varvec{\alpha }}})= \big (\sqrt{2 I_1}\cos (\alpha _1),\,\sqrt{2 I_1}\sin (\alpha _1),\, \sqrt{2 I_2}\cos (\alpha _2),\, \sqrt{2 I_2}\sin (\alpha _2)\big ). \end{aligned}$$

(58)

By applying the Exchange Theorem (see Gröbner 1967; Giorgilli 2003), the solutions of the equations of motions related to $\mathcal {H}_{2DOF}$ can be mapped to those for $\mathcal {H}_{sec, \, 2+3/2}\,$. Indeed, assume that $t\mapsto \big (\tilde{{{\varvec{p}}}}(t),\,\tilde{{{\varvec{q}}}}(t), \,\tilde{{{\varvec{I}}}}(t),\,\tilde{{{\varvec{\alpha }}}}(t)\big )$ is an orbit corresponding to the integrable flow induced by $\mathcal {H}_{2DOF}\,$; in particular, in our model we have that $\tilde{{{\varvec{q}}}}(t)={{\varvec{\omega }}}_Bt={{\varvec{\omega }}}t$, where the components of the angular velocity vector ${{\varvec{\omega }}}$ are given in Eq. (8). Therefore, the orbit

$$\begin{aligned} t\mapsto \big ({{\varvec{\omega }}}t\,, \, {{\mathcal {A}}}\big (\mathscr {F}({{\varvec{\omega }}}t,\tilde{{{\varvec{I}}}}(t),\tilde{{{\varvec{\alpha }}}}(t))\big )\big ) \end{aligned}$$

(59)

is an approximate²² solution of the Hamilton’s Eq. (15).

For our purposes, it is also useful to construct the inverse of the function $\mathscr {F}$, which maps from the original canonical coordinates to the ones referring to the resonant normal form. Therefore, it is convenient to replace all the compositions of Lie series appearing in the r.h.s. of (56) with the following expressions, $\forall \ i=1,2\,$:

$$\begin{aligned} \begin{aligned}&\hat{I}_i=\exp \left( L_{-\chi _0^{(1)}}\right) \exp \left( L_{-\chi _1^{(1)}}\right) \exp \left( L_{-\chi _2^{(1)}}\right) \ldots \\&\quad \ldots \exp \left( L_{-\chi _0^{(\mathcal {N}_S)}}\right) \exp \left( L_{-\chi _1^{(\mathcal {N}_S)}}\right) \exp \left( L_{-\chi _2^{(\mathcal {N}_S)}}\right) I_i\,,\\&\hat{\alpha }_i=\exp \left( L_{-\chi _0^{(1)}}\right) \exp \left( L_{-\chi _1^{(1)}}\right) \exp \left( L_{-\chi _2^{(1)}}\right) \ldots \\&\quad \ldots \exp \left( L_{-\chi _0^{(\mathcal {N}_S)}}\right) \exp \left( L_{-\chi _1^{(\mathcal {N}_S)}}\right) \exp \left( L_{-\chi _2^{(\mathcal {N}_S)}}\right) \alpha _i; \end{aligned} \end{aligned}$$

(60)

gathering all the corresponding changes of coordinates allows us to define²³$\mathscr {C}^{-1}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})$. Proceeding in an analogous way, we can introduce the inverse function of $\mathscr {C}_B\,$; in more detail, we can start from formula (55), by reversing the order of all the Lie series and by changing the sign to all the generating functions, then we can define $\mathscr {C}_B^{-1}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})$. Therefore, we can introduce also

$$\begin{aligned} \mathscr {F}^{-1}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})= \mathscr {C}_B^{-1}\big ({{\varvec{q}}},\mathscr {C}^{-1}({{\varvec{q}}},{{\varvec{I}}},{{\varvec{\alpha }}})\big ). \end{aligned}$$

(61)

We are now ready to exploit the (cheap) numerical solutions of the 2 DOF integrable Hamiltonian, which is described in (53), in order to retrieve information about the secular dynamics of $\upsilon$-And $\textbf{b}$ through our SQPR model. This can be done thanks to the knowledge of the approximate solution (59). The initial conditions are selected in the same way as in Sect. 3.1.2: we consider the initial orbital elements reported in Table 7 and the minimal possible value of the mass of $\upsilon$-And $\textbf{b}\,$, i.e., $m_1=0.674$ $M_{J}\,$. These data are completed with the values of $(\textrm{i}_1(0)\,,\Omega _1(0))$ ranging in the $20\times 60$ regular grid that covers $I_\textrm{i}\times I_\Omega =[6.865^\circ , 34^\circ ]\times [0^\circ ,360^\circ ]$; moreover, all these initial values of the orbital elements are translated in the Laplace reference frame, which refers only to the two outermost exoplanets. Hence, we can compute a set of $21\times 60$ initial conditions of type $\big ({{\varvec{I}}}(0),{{\varvec{\alpha }}}(0)\big )=\mathcal{A}^{-1}\big (\xi _1(0),\eta _1(0),P_1(0),Q_1(0)\big )$, by using formula (2) with $j=1\,$, (50) and the definition (58). Finally, we can translate the initial conditions to initial values of the canonical coordinates found after the resonant normal form, by computing $\big (\tilde{{{\varvec{I}}}}(0),\tilde{{{\varvec{\alpha }}}}(0)\big )= \mathscr {F}^{-1}\big ({{\varvec{0}}},{{\varvec{I}}}(0),{{\varvec{\alpha }}}(0)\big )$.

As shown below, an important information is obtained by a criterion allowing to identify those domains of initial conditions in which the series are either divergent or slowly converging. We introduce such a criterion to preselect initial conditions that are admissible. From a mathematical point of view, the identity $({{\varvec{I}}}(0),{{\varvec{\alpha }}}(0))=({{\varvec{I}}}^{O}(0),{{\varvec{\alpha }}}^{O}(0)):= \mathscr {F}\big ({{\varvec{0}}},\mathscr {F}^{-1}\big ({{\varvec{0}}},{{\varvec{I}}}(0),{{\varvec{\alpha }}}(0)\big )\big )\,$ holds in a domain where the normalization procedure is convergent, provided that no truncations are applied to the series $\mathscr {F}$ and $\mathscr {F}^{-1}$ and that the computation of the series is not affected by any round-off errors. Due to the errors and truncations introduced in the computation, however, in general we obtain that $({{\varvec{I}}}(0),{{\varvec{\alpha }}}(0))\ne ({{\varvec{I}}}^{O}(0),{{\varvec{\alpha }}}^{O}(0))\,$. In the domain where the series expansions are rapidly converging the difference $({{\varvec{I}}}(0),{{\varvec{\alpha }}}(0))-({{\varvec{I}}}^{O}(0),{{\varvec{\alpha }}}^{O}(0))$ is small. When, instead, we obtain a large difference, this is an indicator that we are outside the domain of convergence of the series. The situation is represented graphically in Fig. 7.

Fig. 7

Graphical representation of the definitions about the initial conditions

In view of the above, we define the following preselection criterion of admissible initial conditions. For any initial condition $({{\varvec{I}}}(0),{{\varvec{\alpha }}}(0))$, we compute the quantities

$$\begin{aligned} \mathfrak {r}_1=\frac{\sqrt{I_1(0)}-\sqrt{I_1^{O}(0)}}{\sqrt{I_1(0)}}\,, \qquad \quad \mathfrak {r}_2=\frac{\sqrt{I_2(0)}-\sqrt{I_2^{O}(0)}}{\sqrt{I_2(0)}}\,. \end{aligned}$$

(62)

The use of the quantities $\sqrt{I_1}$ and $\sqrt{I_2}$ is motivated by the fact that they are of the same order of magnitude as the eccentricity and the inclination of $\upsilon$-And $\textbf{b}\,$, respectively. Moreover, it is useful to define also the following ratios

$$\begin{aligned} \mathfrak {R}_1(t)= \frac{\sqrt{\big (\tilde{\xi }_1(t)\big )^2+\big (\tilde{\eta }_1(t)\big )^2}}{\sqrt{\big (\tilde{\xi }_1(0)\big )^2+\big (\tilde{\eta }_1(0)\big )^2}} \,,\qquad \qquad \mathfrak {R}_2(t)= \frac{\sqrt{\big (\tilde{P}_1(t)\big )^2+\big (\tilde{Q}_1(t)\big )^2}}{\sqrt{\big (\tilde{P}_1(0)\big )^2+\big (\tilde{Q}_1(0)\big )^2}}\,, \end{aligned}$$

(63)

where

$$\begin{aligned} t\mapsto \big ({{\varvec{\omega }}}t,\,\tilde{\xi }_1(t),\tilde{\eta }_1(t),\tilde{P}_1(t),\tilde{Q}_1(t)\big ):= \big ({{\varvec{\omega }}}t,\, \mathcal{A}\big (\mathscr {F}({{\varvec{\omega }}}t,\tilde{{{\varvec{I}}}}(t), \tilde{{{\varvec{\alpha }}}}(t))\big )\big ) \end{aligned}$$

is the approximate solution of Hamilton’s equations (15), as produced by the semi-analytic integration scheme summarized in formula (59). Comparing formula (63) with the definition of the Poincaré canonical variables in (2), it is easy to realize that $\mathfrak {R}_1$ and $\mathfrak {R}_2$ are functions of the time that describe the behavior of the orbital excursions with respect to the eccentricity and the inclination of $\upsilon$-And $\textbf{b}\,$, respectively. We then investigate the behavior of the following function:

$$\begin{aligned} \tilde{\textrm{e}}_1(t)=\sqrt{\frac{2 \tilde{I}_1(t)}{\Lambda _1} - \frac{\tilde{I}_1^2(t)}{\Lambda _1^2}}\,. \end{aligned}$$

(64)

Note that $\tilde{\textrm{e}}_1$ would be equal to the eccentricity of $\upsilon$-And $\textbf{b}$ if $\tilde{I}_1$ was replaced by $\big (\xi _1^2+\eta _1^2\big )/2$, with $(\xi _1,\,\eta _1)$ defined in (2). However, the new action $\tilde{I}_1$ is conjugated to $\big (\tilde{\xi }_1^2+\tilde{\eta }_1^2\big )/2$ which is only nearly equal to $\big ({\xi }_1^2+{\eta }_1^2\big )/2\,$, since the composition of the transformations $\mathscr {C}$ and $\mathscr {C}_B$ is near-to-identity. Therefore, we can consider $\tilde{\textrm{e}}_1$ as an approximate evaluation of the eccentricity under the resonant normal form model.

For each pair $\big (\textrm{i}_1(0),\Omega _1(0)\big )\,$ of the $21\times 60$ points defining the grid which covers $I_\textrm{i}\times I_\Omega =[6.865^\circ , 34^\circ ]\times [0^\circ ,360^\circ ]$ we determine the corresponding initial conditions of type $({{\varvec{I}}}(0),{{\varvec{\alpha }}}(0))$, as explained above, and we proceed as follows:

• if $\max \{\mathfrak {r}_1,\,\mathfrak {r}_2\}> 1$, then the corresponding initial condition is considered as “non-admissible,” i.e., outside the domain of applicability of the series. Then, we skip the step below and pass directly to consider the next initial conditions of the grid;

• If the initial conditions is admissible, we numerically²⁴ solve the equations of motion for the integrable Hamiltonian model with 2 DOF described in (53), using a RK4 method and starting from $\big ({{\varvec{0}}},\,\tilde{{{\varvec{I}}}}(0),\tilde{{{\varvec{\alpha }}}}(0)\big )$; during such a numerical integration, we compute the maximal values attained by the three previously defined quantitative indicators, that are

  $$\begin{aligned} \mathfrak {R}_{1\,\textrm{MAX}} = \max _t\{\mathfrak {R}_1(t)\}, \quad \mathfrak {R}_{2\,\textrm{MAX}} = \max _t\{\mathfrak {R}_2(t)\}, \quad \tilde{\textrm{e}}_{1\,\textrm{MAX}} = \max _t\{\tilde{\textrm{e}}_1(t)\}. \end{aligned}$$

Fig. 8

Color-grid plots of the maximal values reached by the ratio $\mathfrak {R}_1$ (on the left) and $\mathfrak {R}_2\,$ (on the right); see the text for more details

Fig. 9

Color-grid plots of the maximum of the function $\tilde{\textrm{e}}_{1}(t)$, which is defined in (64)

The results about the maxima of the functions defined in (63)–(64) are reported in Figs. 8 and 9. The white central regions of those pictures correspond to those pairs $\big (\textrm{i}_1(0),\Omega _1(0)\big )\,$ for which we obtain failure of the preliminary test, i.e., $\max \{\mathfrak {r}_1,\,\mathfrak {r}_2\}>1$. We immediately recognize that the missing part of the plots (where the determination of the initial conditions is considered so unreliable that the corresponding numerical integrations are not performed at all) nearly coincides with the central region of Fig. 2a, where the orbital eccentricity of $\upsilon$-And $\textbf{b}$ reaches critical values. We conclude that the stability domain in the space of the initial values of $\textrm{i}_1(0)$ and $\Omega _1(0)$ (which are unknown observational data) can be reconstructed in a reliable way through the application of the above criterion, which only involves the series transformations, as well as through the numerical solutions of our integrable secular model with 2 DOF. We emphasize that this allows to reduce significantly the computational cost with respect to the long-term symplectic integrations of the complete four-body problem, which is a 9 DOF Hamiltonian system.

Comparing the regions of the stability domain at the border near the (white) central ones, we see that all three numerical indicators plotted in Figs. 8 and 9 increase their values when the unstable zone is approached. This is in agreement with the expectations and the comparison with Fig. 2a. On the other hand, the 2 DOF secular model is unable to capture the region of instability internal to the stable one, highlighted by two green stripes starting from the bottom of Fig. 2a in correspondence with $\Omega _1(0)=0^\circ =360^\circ$. The two curved stripes look rather symmetric, and they join each other around the point $\big (\textrm{i}_1(0),\,\Omega _1(0)\big )\simeq \big (30^\circ ,\,0^\circ \big )=\big (30^\circ ,\,360^\circ \big )$. Since the dependence of the Hamiltonian on the angles ${{\varvec{q}}}$ (which describes the dynamics of the outer exoplanets) is removed from the 2DOF model, it seems reasonable that some of the resonances are not present in the normal form generated by the algorithm à la Birkhoff, although they play a remarkable role in the dynamics of more complex models.

6 Secular orbital evolution of $\upsilon$-And $\textbf{b}$ taking also into account relativistic effects

In this section, we study the dynamics of $\upsilon$-And $\textbf{b}$ in the framework of a secular quasi-periodic restricted Hamiltonian model where also corrections due to general relativity are taken into account. Since we focus on the orbital dynamics of the innermost planet of the $\upsilon$-Andromedæ system and it is very close to a star that is about 30% more massive than the Sun (let us recall that the value of the semi-major axis of $\upsilon$-And $\textbf{b}$ is reported in Table 7, i.e., $a_1=0.0594$ AU), it is natural to expect that the corrections due to general relativity can play a relevant role for the system under consideration. Similarly as in the previous sections, we study these effects in the framework of a 2 DOF secular model. We start by considering the following Hamiltonian:

$$\begin{aligned} \mathcal {H}=\mathcal {H}_{4BP}+\mathcal {H}_{GR}\,, \end{aligned}$$

where $\mathcal {H}_{4BP}$ defines the four-body problem (see (9)) and $\mathcal {H}_{GR}$ describes the general (post-Newtonian) relativistic corrections to the Newtonian mechanics. Following Migaszewski and Goździewski (2009), the secular quasi-periodic restricted Hamiltonian which includes corrections due to the general relativity (hereafter, GR) is obtained by removing the dependence of the Hamiltonian on the fast angles. Therefore, we introduce

$$\begin{aligned} \mathcal {H}^{(GR)}_{sec}= \int _{\mathbb {T}^3}\frac{\mathcal {H}_{4BP}}{8\pi ^3}\, \textrm{d}\lambda _1 \textrm{d}\lambda _2 \textrm{d}\lambda _3 + \int _{\mathbb {T}} \frac{\mathcal {H}_{GR}}{2\pi }\, \textrm{d} M_1\,:= \mathcal {H}^{(NG)}_{sec}+\left\langle \mathcal {H}_{{GR}} \right\rangle _{M_1} \,, \end{aligned}$$

(65)

where the expansion of the mean of the 4BP Hamiltonian $\mathcal {H}^{(NG)}_{sec}$ (recall definition (10)) is explicitly written in Eq. (11), while the average of the GR contribution with respect to the mean anomaly of $\upsilon$-And $\textbf{b}$ is given by

$$\begin{aligned} \left\langle \mathcal {H}_{{GR}} \right\rangle _{M_1}= -\frac{3\,\mathcal {G}^2\, m_0^2\, m_1}{ a_1^2 c^2\sqrt{1-\textrm{e}_1^2}} +\frac{15\,\mathcal {G}^2\, m_0^2\, m_1}{8 a_1^2 c^2} -\frac{\mathcal {G}^2 \,m_0\, m_1^2}{8 a_1^2 c^2}\,, \end{aligned}$$

(66)

c being the velocity of light in vacuum. In the above expression of $\left\langle \mathcal {H}_{{GR}} \right\rangle _{M_1}\,$, the summand where the eccentricity of $\upsilon$-And $\textbf{b}$ (i.e., $\textrm{e}_1$) occurs in the denominator is the only to be untrivial, in the sense that the other two give additional constant contribution to the secular Hamiltonian and, then, they can be disregarded. By proceeding in a similar way to what has been already done for the classical expansions of the initial Hamiltonian (1), it is possible to express $\left\langle \mathcal {H}_{{GR}} \right\rangle _{M_1}$ in the Poincaré variables $(\xi _1,\eta _1)$, described in Eq. (2).

Thus, keeping in mind the procedure explained in Sect. 3, one easily realizes that the secular quasi-periodic restricted model of the dynamics of $\upsilon$-And $\textbf{b}$ which includes relativistic corrections (hereafter, SQPR-GR) can be described by the following $2+3/2$ DOF Hamiltonian:

$$\begin{aligned} \begin{aligned} \mathcal {H}^{(GR)}_{sec,\, 2+\frac{3}{2}}&({{\varvec{p}}},{{\varvec{q}}}, \xi _1, \eta _1, P_1, Q_1 )= \omega _3\,p_3 + \omega _4\,p_4 + \omega _5\,p_5\\&\qquad +\mathcal {H}^{(NG)}_{sec}(q_3, q_4, q_5, \xi _1, \eta _1, P_1, Q_1 ) +\left\langle \mathcal {H}_{GR} \right\rangle _{M_1}(\xi _1,\eta _1) \,, \end{aligned} \end{aligned}$$

(67)

where the angular velocity vector ${{\varvec{\omega }}}=(\omega _3,\omega _4,\omega _5\,)$ is given in (8) and $\mathcal {H}^{(NG)}_{sec}$ can be replaced by $\mathcal {H}^{1-2}_{sec}+\mathcal {H}^{1-3}_{sec}$ appearing in formula (13). Finally, in the framework of this SQPR-GR model, the equations for the orbital motion of the innermost planet can be written as

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q_3}=\partial \mathcal {H}^{(GR)}_{sec,2+\frac{3}{2}}/\partial p_3=\omega _3\\ \dot{q_4}=\partial \mathcal {H}^{(GR)}_{sec,2+\frac{3}{2}}/\partial p_4=\omega _4\\ \dot{q_5}=\partial \mathcal {H}^{(GR)}_{sec,2+\frac{3}{2}}/\partial p_5=\omega _5\\ \dot{\xi }_{1}= -\partial \left( \mathcal {H}_{sec}^{(NG)}(q_3, q_4, q_5, \xi _1, \eta _1, P_1, Q_1 ) +\left\langle \mathcal {H}_{GR} \right\rangle _{M_1}(\xi _1, \eta _1 )\right) /\partial \eta _1 \\ \dot{\eta }_{1}= \partial \left( \mathcal {H}_{sec}^{(NG)}(q_3, q_4, q_5, \xi _1, \eta _1, P_1, Q_1 ) +\left\langle \mathcal {H}_{GR} \right\rangle _{M_1}(\xi _1, \eta _1)\right) /\partial \xi _1 \\ \dot{P}_{1}= -\partial \mathcal {H}_{sec}^{(NG)}(q_3, q_4, q_5, \xi _1, \eta _1, P_1, Q_1) /\partial Q_1 \\ \dot{Q}_{1}= \partial \mathcal {H}_{sec}^{(NG)}(q_3, q_4, q_5, \xi _1, \eta _1, P_1, Q_1 ) /\partial P_1 \end{array}\right. }\,. \end{aligned}$$

(68)

6.1 Numerical integration of the SQPR-GR model

Similarly as in Sect. 3.1.2, we numerically integrate the equations of motion for the secular quasi-periodic restricted Hamiltonian with general relativistic corrections, defined in formula (68). As initial values of the orbital parameters $a_1(0)$, $\textrm{e}_1(0)$, $M_1(0)$ and $\omega _1(0)$ we take those reported in Table 7; moreover, we set $m_1=0.674$ as value for the mass of $\upsilon$-And $\textbf{b}$ and $(\textrm{i}_1(0)\,,\Omega _1(0))$ ranging in the $20\times 60$ regular grid that covers $I_\textrm{i}\times I_\Omega =[6.865^\circ , 34^\circ ]\times [0^\circ ,360^\circ ]$. Hence, it is possible to compute the corresponding initial values of the orbital elements in the Laplace reference frame (which is determined taking into account $\upsilon$-And $\textbf{c}$ and $\upsilon$-And $\textbf{d}$ only) and to perform $21\times 60$ numerical integrations starting from all these initial data. Once again, for each numerical integration, we are interested in determining the maximal values reached by the eccentricity of $\upsilon$-And $\textbf{b}$ and by the maximal mutual inclination between $\upsilon$-And $\textbf{b}$ and $\upsilon$-And $\textbf{c}\,$. The results are reported in the color-grid plots of Fig. 10.

Fig. 10

Color-grid plots of the maximal value reached by the eccentricity of $\upsilon$-And $\textbf{b}$ (on the left) and by the mutual inclination between $\upsilon$-And $\textbf{b}$ and $\upsilon$-And $\textbf{c}$ (on the right). The maxima are computed during the RK4 numerical integrations (each of them covering a timespan of $10^5$ yr) of the SQPR-GR equations of motion (68), which cover a timespan of $10^5$ yr

Fig. 11

Behavior of the fundamental angular velocity $\nu$ as obtained by applying the frequency map analysis method to the signal $\xi _1(t)+i\eta _1(t)\,$, computed through the RK4 numerical integration of the SQPR-GR model (68), covering a timespan of $1.31072\cdot 10^{5}$ yr. We take, as initial conditions, $(\textrm{i}_1(0),\Omega _1(0))\in \lbrace 6.865^\circ \rbrace \times I_\Omega$ for the left panel and $(\textrm{i}_1(0),\Omega _1(0))\in \lbrace 10.9353^\circ \rbrace \times I_\Omega$ for the right one

By comparing Fig. 10a with Fig. 5a, one can immediately realize that the regions colored in blue are much wider in the former than in the latter one. Indeed, the darker regions refer to initial conditions which generate motions with maximal values of the eccentricity of $\upsilon$-And $\textbf{b}$ that are relatively low, while the zones colored in red or yellow correspond to such large values of the eccentricity implying that those orbits have to be considered unstable. Therefore, our numerical explorations highlight that the effects due to general relativity play a stabilizing role on the orbital dynamics of the innermost planet. This conclusion is in agreement with was already remarked about the past evolution of our Solar System, in particular for what concerns the orbital eccentricity of Mercury (see Laskar and Gastineau 2009).

Moreover, as already done in Sect. 3.1.2, in order to further explore the stable and chaotic regions of Fig. 10a, we apply the frequency map analysis method to the signal $\xi _1(t)+i\eta _1(t)$ as produced by the numerical integration of the system (68), i.e., in the SQPR-GR approximation. We perform the numerical integrations as described at the beginning of the present Section, taking into account only a few values in $I_\textrm{i}\,$ for the initial inclinations, i.e., $\textrm{i}_1(0)=6.865^\circ ,\, 8.22175^\circ ,\,9.5785^\circ ,\,10.9353^\circ$ and $\Omega _1(0)\in I_\Omega \,$. In Fig. 11 we report the behavior of the angular velocity $\nu$ corresponding to the first component of the approximation of $\xi _1(t)+i\eta _1(t)$, as obtained by applying the FA computational algorithm; we recall that this quantity is related to the precession rate of $\varpi _1\,$. As initial value for the inclination $\textrm{i}_1(0)$ we fix $6.865^\circ$ for Fig. 11a and $10.9353^\circ$ for Fig. 11b. Also here, we do not report the cases $(\textrm{i}_1(0),\Omega _1(0))\in \lbrace 8.22175^\circ ,\, 9.5785^\circ \rbrace \times I_\Omega \,$, since the behavior of these plots is similar to the ones in Fig. 11.

The situation is well described by Fig. 11a and analogous considerations hold for Fig. 11b apart a few main differences which will be highlighted in the following discussion. When the values of $\Omega _1(0)$ are ranging in $[0, \sim 120^\circ ]$ and $[\sim 260^\circ , \, 360^\circ ]$ we can observe a regular behavior of the angular velocity $\nu$ which is also nearly monotone, with the only exception around a local minimum. According to the frequency map analysis method, such a regular regime is due to the presence of many invariant tori which fill the stability region located at the two lateral sides of the plot 10a. In the case of Fig. 11a, this also applies when $\Omega _1(0)$ is ranging in $[\sim 150^\circ , \sim 220^\circ ]$, which corresponds to the stable blue internal area of Fig. 10a. On the other hand, in the case of Fig. 11b, for the same range of initial values of the node longitude of $\upsilon$-And $\textbf{b}$, the behavior is not so regular; this is in agreement with the fact that in correspondence with $\textrm{i}_1(0)\sim 11^\circ$ the plot of the maximal values of $\textrm{e}_1$ in the central region highlights the occurrence of chaotical phenomena. Moreover, for what concerns values of $\Omega _1(0)$ in $[\sim 120^\circ , \, \sim 150^\circ ]$ and $[\sim 220^\circ , \, \sim 260^\circ ]$ (corresponding to the green stripes of Fig. 10a), Fig. 11a shows a behavior typical of the crossing of a resonance in the chaotic region surrounding a separatrix. The value of the angular velocity for which this phenomenon takes place is, again, related to $\omega _4\simeq -1.04\cdot 10^{-3}$ (as it can be easily appreciated looking to the small plateau appearing in Fig. 11b).

Comparing Fig. 11 with Fig. 6, the enlargement of the stable region is evident. Moreover, we can also see how much this phenomenon is influenced by the modification of the pericenter precession rate of the inner planet due to relativistic effects. Indeed, looking at the values reported on the y-axis of Figs. 11 and 6, one can appreciate that the fundamental angular velocity, in the case of the SQPR model, takes values remarkably closer to zero with respect to those assumed in the case of the SQPR-GR model.

6.2 Application of the normalization algorithms to the secular quasi-periodic restricted model of the dynamics of $\upsilon$-And $\textbf{b}$ with relativistic corrections

Starting from Hamiltonian (67), we can reapply the normalization algorithms described in Sects. 4.1 and 4.2. All this computational procedure ends up with the introduction of a new 2 DOF Hamiltonian²⁵ model which can be written in the following form (analogous to the one reported in formula (53)):

$$\begin{aligned} \begin{aligned} \mathcal {H}_{2DOF}^{(GR)}({{\varvec{I}}},{{\varvec{\alpha }}})&= \mathcal {E}_{B;GR}+{{\varvec{\Omega }}}_{B;GR}\cdot {{\varvec{I}}}+ \sum _{s= 0}^{\mathcal {N}_S}\sum _{l=3}^{\mathcal {N}_L} h_{l}^{(\mathcal {N}_L-2,s)}({{\varvec{I}}},{{\varvec{\alpha }}})\, , \end{aligned} \end{aligned}$$

(69)

where $\mathcal {E}_{B;GR}\in \mathbb {R}$, ${{\varvec{\Omega }}}_{B;GR}\in \mathbb {R}^2$ and $h_{l}^{(\mathcal {N}_L-2,s)}\in \mathfrak {P}_{l,2\,s}$ $\forall \ l=3,\,\ldots \,,\,\mathcal {N}_L\,,\> s=0\,,\,\ldots \,,\,\mathcal {N}_S\,$. We emphasize that also $\mathcal {H}_{2DOF}^{(GR)}$ is integrable because of the same reasons already discussed in Sect. 5; indeed, after having checked that ${\partial \mathcal {H}_{2DOF}^{(GR)}}/{\partial \alpha _1}\,+\, {\partial \mathcal {H}_{2DOF}^{(GR)}}/{\partial \alpha _2}\,=\,0\,$, we can apply the Liouville theorem, because there are two independent constants of motion, i.e., $I_1+I_2$ and $\mathcal {H}_{2DOF}^{(GR)}$ itself.

Moreover, also for this new model we can reproduce the same kind of numerical exploration described in Sect. 5. In particular, we can compute the values of the numerical indicators $\mathfrak {R}_{1\,\textrm{MAX}}\,$, $\mathfrak {R}_{2\,\textrm{MAX}}$ and $\tilde{\textrm{e}}_{1\,\textrm{MAX}}$ corresponding to each pair $\big (\textrm{i}_1(0),\Omega _1(0)\big )\,$ of the $21\times 60$ points defining the regular grid which covers $I_\textrm{i}\times I_\Omega =[6.865^\circ , 34^\circ ]\times [0^\circ ,360^\circ ]$. In the following, we analyze the color-grid plots for a few different values of the parameter ruling the truncation in the trigonometric degree, namely $\mathcal {N}_S\,$, and in the square root of the action, i.e., $\mathcal {N}_L\,$. The color-grid plots for the maximal value reached by $\mathfrak {R}_1\,$ and $\tilde{\textrm{e}}_1$ are reported in Figs. 12, 13 and 14.

Fig. 12

Color-grid plots of the maximal values reached by the ratio $\mathfrak {R}_1(t)$ (on the left) and the function $\tilde{\textrm{e}}(t)$ (on the right), which are defined in (63)–(64). These laws of motion are computed along the flow induced by the 2 DOF Hamiltonian $\mathcal {H}_{2DOF}^{(GR)}$ which takes into account also GR corrections and is defined in (69), in the particular case with $\mathcal {N}_S=5$ and $\mathcal {N}_L=6$

Fig. 13

Same as in Fig. 12, in the case with $\mathcal {N}_S=6$ and $\mathcal {N}_L=5$

Let us recall that $\mathfrak {R}_{2\,\textrm{MAX}}$ is an evaluation of the maximal excursion of the inclination of $\upsilon$-And $\textbf{b}$. For the sake of brevity, its plots are omitted and they are not included in Figs. 12, 13 and 14, because in our numerical explorations the ranges of values experienced by $\mathfrak {R}_{2\,\textrm{MAX}}$ are so narrow that their analysis does not look so significant. Therefore, it is better to focus on the plots of $\mathfrak {R}_{1\,\textrm{MAX}}$ and $\tilde{\textrm{e}}_{1\,\textrm{MAX}}\,$; let us recall that both of them refer to the eccentricity of $\upsilon$-And $\textbf{b}$. By comparing Figs. 12, 13 and 14, one can appreciate a well-known phenomenon concerning the constructive algorithms à la Birkhoff: the greater the number of normalization steps (i.e., $\mathcal {N}_L-2$), the smaller the domain of applicability (see, e.g., Giorgilli 2003 for the discussion about the determination of the optimal step).

Fig. 14

Same as in Fig. 12, in the case with $\mathcal {N}_S=6$ and $\mathcal {N}_L=4$

By comparing Figs. 13 and 14 also with Fig. 10a, we observe that in the cases with $\mathcal {N}_L=4,\,5$ our computational procedure is able to reconstruct with a good accuracy the U-shaped border of the stability domain. Note that the horizontal strip at the bottom of the plots²⁶ corresponds to orbital motions which look stable, since the eccentricity of $\upsilon$-And $\textbf{b}$ does not reach large values (with the eventual exception of the narrow green areas that in Fig. 10a are expected to correspond to resonant regions). This highlights a main difference with the non relativistic model discussed in Sect. 5, because in that case there is an interval of values of $\Omega _1(0)$ centered about $180^\circ$ for which none of the initial inclinations $\textrm{i}_1(0)\in [6.865^\circ , 34^\circ ]$ corresponds to a stable orbital configuration (see Fig. 5a). The reliability of our simplified 2 DOF Hamiltonian model (which is defined in formula (69) and takes into account also GR corrections) is enforced also by the fact that it is able to capture also this phenomenon.